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Abstract 

This  research  modeled  low-speed  flow  past  idealized  engine  nacelle  clutter  in  sup¬ 
port  of  aircraft  fire  suppression  research.  The  idealized  clutter  was  comprised  of  three 
vertical  rows  of  staggered  circular  cylinders  approximating  typical  nacelle  obstruc¬ 
tions  such  as  fuel  lines  and  wire  bundles.  Single-phase,  Detached-Eddy  Simulations 
(DES)  were  conducted  using  the  commercial  CFD  solver,  Fluent^ ,  to  resolve  the 
flow-field  dynamics  inside  the  clutter  element  and  determine  mechanisms  accounting 
for  the  failure  of  suppressant  spray  droplets  from  traversing  the  array  under  low-speed, 
free-stream  conditions  (Rep  =  1,575).  The  numerical  models  provided  no  evidence 
that  span-wise  vorticity  or  non-uniform  shedding  was  responsible  for  transporting 
dispersed-phase  particles  towards  the  tunnel  walls  for  deposition.  However,  the  simu¬ 
lations  demonstrated  that  suppressant  droplets  would  likely  follow  a  path  governed  by 
the  vector  sum  of  the  local  carrier  fluid  velocity  and  the  velocity  imposed  by  gravity. 
Additionally,  the  Stokes  number  was  computed  from  time-accurate  data  to  determine 
the  ability  of  dispersed  particles  to  negotiate  the  clutter  element  without  impinging 
on  a  cylinder.  For  slower  free-stream  velocities,  U0 0  =  1  m/s,  suppressant  droplets 
(D  =  90  /mi)  will  likely  be  entrained  in  vortices  shed  from  the  intermediate  row 
of  cylinders  and  subsequently  deposited  on  the  last  row  of  cylinders  as  the  Karman 
vortex  directly  collides  with  the  clutter.  At  free-stream  velocities,  =  5  m/s,  the 
droplet  particles  will  likely  fail  to  track  the  carrier  fluid  streamlines  in  the  cylinder 
wake  and  remain  free  of  any  shed  vortices.  Thus,  the  suppressant  will  conceivably 
transit  the  cylinder  array  without  impact.  These  findings  imply  that  a  bluff-body 
turbulent  diffusion  flame  in  a  cylinder  wake  could  be  nearly  impossible  to  extinguish 
under  high-speed,  co-flow  conditions.  Conversely,  suppressant  transported  by  low- 
speed  co-flow  would  experience  difficulty  traversing  the  cylinder  array  and  reaching  a 
downstream  fire. 
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CFD  Investigation  of  Flow  Past 


Idealized  Engine  Nacelle  Clutter 

I.  Introduction 

1.1  Engine  Nacelle  Fire  Suppression 

The  aircraft  engine  nacelle  environment  is  particularly  conducive  to  the  outbreak 
of  fire  as  a  result  of  the  myriad  of  ignition  sources  such  as  chaffed  electrical 
wiring,  hot  surfaces,  and  the  abundance  of  fuel  from  leaking  or  damaged  jet  fuel,  oil, 
or  hydraulic  lines.  Ambient  air  typically  enters  the  engine  nacelle  through  vents  or 
actuated  doors  and  is  used  to  cool  the  engine  core  and  avionics,  as  well  as  disperse 
flammable  vapors  which  would  otherwise  linger  near  hot  surfaces.  In  the  event  of  a 
fire,  the  flow  of  cooling  air  through  the  nacelle  has  the  undesirable  effect  of  diluting 
a  dispensed  fire  suppressant  and  raising  the  probability  of  re-ignition.  As  cooling  air 
flows  through  the  nacelle  it  is  forced  to  navigate  a  labyrinth  of  clutter  in  the  form  of 
wire  bundles,  fuel  lines,  oil  pumps,  avionics  boxes,  and  structural  members  such  as 
ribs  and  stringers.  The  obstructions  give  rise  to  all  types  of  flow  irregularities  in  an 
already  highly  turbulent  free-stream  flow  field.  Given  these  conditions,  the  most  likely 
type  of  nacelle  fire  is  either  a  turbulent  diffusion  flame  stabilized  by  a  clutter  element 
or  a  pool  fire  guarded  by  an  obstruction.  [1]  Naturally,  the  ability  to  extinguish  a 
nacelle  fire  is  paramount  to  the  safety  and  survivability  of  the  aircraft  and  crew. 

In  use  since  the  1950s,  Halon  1301  (CF3Br)  is  operationally  proven  to  be  the 
fire  suppressant  of  choice  for  both  military  and  commercial  aircraft  engine  nacelles. 
As  a  result  of  the  Montreal  Protocol  on  Substances  that  Deplete  the  Ozone  Layer, 
the  production  of  Halon  was  banned  as  of  January  1,  1994  forcing  the  Department 
of  Defense  to  rely  on  its  stockpiles  of  Halon  1301  for  aircraft  fire  suppression.  [13] 
The  ban  sparked  the  creation  of  the  LLS.  Air  Force  Halon  Replacement  Program  for 
Aviation  and  the  Department  of  Defense  Next  Generation  Fire  Suppression  Technol- 
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ogy  Program  (NGP);  both  programs  were  responsible  for  evaluating  alternative  ex- 
tinguishants  and  developing  technology  demonstrations  of  environmentally-friendly, 
economical  methods  for  fire  suppression  and  halon  replacement.  The  NGP  primary 
areas  of  study  included  fire  extinguishant  chemistry  and  toxicity  testing,  developing 
new  aerosol  suppressants,  improving  methods  of  suppressant  delivery,  and  the  ex¬ 
amination  of  new  suppression  mechanisms  and  technology.  [17]  Because  Halon  1301 
extinguishes  aircraft  fires  so  efficiently,  research  to  characterize  the  flow  held  and  opti¬ 
mize  suppressant  delivery  had  little  pay-off  until  the  production  ban  and  the  initiation 
of  the  NGP.  [7] 

At  the  outset  of  the  NGP  and  the  Air  Force  Halon  Replacement  Program  for 
Aviation,  extensive  testing  was  conducted  at  the  Aircraft  Engine  Nacelle  Fire  Test 
Simulator  (AENFTS)  facility  located  at  Wright- Patterson  Air  Force  Base  in  Dayton, 
Ohio.  The  AENFTS  is  constructed  to  represent  the  geometry  and  how  helds  associ¬ 
ated  with  aircraft  engine  nacelles.  As  part  of  the  NGP,  factors  such  as  extinguishant 
temperature,  pressure,  and  distribution,  ventilation  air  pressure  and  velocity,  clutter 
size  and  conhguration,  and  hre  location  were  examined  at  the  AENFTS  facility  and 
correlated  in  an  effort  to  develop  metrics  for  evaluating  hre  suppression.  The  ex¬ 
perimental  characterization  of  nacelle  hres  provided  insight  into  airhow  velocities  and 
initial  estimates  of  minimum  agent  concentrations  necessary  for  hre  suppression.  [1,13] 
However,  the  details  of  agent  transport  were  still  largely  unexamined. 

Working  under  the  suppressant  transport  branch  of  the  NGP,  Disimile  et  al.  [7] 
experimentally  explored  water  droplet  transport  past  idealized  engine  nacelle  clutter. 
Idealized  clutter  models  are  simplified  representations  of  actual  engine  nacelle  clutter 
elements  which  approximate  complicated  nacelle  obstructions  allowing  investigators 
to  isolate  geometric  variables.  For  example,  structural  ribs  are  approximated  as  simple 
fences,  wire  bundles  and  fuel  lines  are  simplified  as  either  a  single  cylinder  or  an  array 
of  cylinders,  and  extremely  dense  areas  of  clutter  are  approximated  as  tightly  packed 
spheres  or  porous  media  (Figure  1.1).  [10] 
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Figure  1.1:  Idealized  Engine  Nacelle  Clutter  Elements  [17] 
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Figure  1.2:  Idealized  Engine  Nacelle  Clutter  Element  of  Disimile  et  al.  [7] 


The  generic  clutter  package  of  Disimile  et  al.  shown  in  Figure  1.2  was  comprised 
of  circular  cylinders  each  50.8  mm  in  diameter  arranged  in  an  array  of  three  equally 
spaced  rows  of  five,  six,  and  five  cylinders.  The  study  examined  the  effect  of  two 
variables,  free-stream  airspeed  and  cylinder  spacing,  on  the  amount  of  suppressant 
transported  past  the  clutter  package.  The  cylinder  spacing  was  varied  from  0.25D  to 
2.0 D,  and  the  free-stream  airspeed  was  examined  over  a  range  of  0.5  mj s  to  6.5  m/s. 
Disimile  et  al.  found  that  for  large  clutter  spacing  (2.0 D)  and  higher  airspeeds  (5  m/s) 
upwards  of  95%  of  the  water  spray  was  transported  past  the  clutter  package.  However, 
for  tight  stream-wise  cylinder  spacing  (0.25D)  and  lower  airspeeds  (0.5—2  m/s)  less 
than  40%  of  the  suppressant  was  transported  past  the  clutter.  [7] 
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1.2  Research  Objective 

Because  of  the  growing  need  to  accurately  understand  flow  field  dynamics  inside 
an  engine  nacelle  for  fire  suppression  applications,  this  research  seeks  to  characterize 
the  unsteady  flow  dynamics  inside  the  idealized  engine  nacelle  clutter  of  Disimile 
et  al.  and  ascertain  the  mechanisms  of  droplet  transport  and  entrapment  observed 
by  Disimile  et  ah  Naturally,  a  deeper  understanding  of  the  processes  contributing 
to  transport  and  collection  of  entrained  droplets  inside  engine  nacelles  is  crucial  to 
formulating  a  comprehensive  nacelle  fire  suppression  strategy. 

1.3  Research  Approach 

Investigation  of  flow  dynamics  inside  the  idealized  clutter  element  of  Disimile  et 
ah  and  the  analysis  of  droplet  transport  phenomena  was  accomplished  by  performing 
numerical  simulations  of  the  flow-held  using  the  commercial  Computational  Fluid  Dy¬ 
namics  (CFD)  solver,  Fluent™.  Detached-Eddy  Simulation  (DES),  a  relatively  recent 
addition  to  the  turbulence  modeling  arsenal,  was  selected  as  the  primary  turbulence 
modeling  technique  because  the  turbulent  structures  most  likely  coupled  to  fire  sup¬ 
pression  and  droplet  transport  are  not  adequately  captured  by  mean  how  quantities. 
Simulations  to  validate  DES  as  an  acceptable  model  in  low  Reynolds  number  hows 
were  conducted  hrst  on  a  single  cylinder  grid.  Subsequently,  DES  was  applied  to 
a  scaled  version  of  the  idealized  nacelle  clutter  developed  by  Disimile  et  ah  at  two 
distinct  free-stream  airspeeds.  Unlike  the  experiment  of  Disimile  et  ah,  the  numerical 
simulation  did  not  include  a  dispersed-phase  model  for  particle  tracking.  Instead,  this 
research  focused  on  resolving  the  how-held  dynamics  and  pinpointing  aspects  of  the 
how  which  might  explain  the  droplet  transport  phenomena. 

Two  hypotheses  explaining  the  transport  and  the  entrapment  of  entrained  sup¬ 
pressant  droplet  particles  were  specifically  tested  in  this  investigation.  One  hypothesis 
postulated  that  non-uniform,  span-wise  shedding  was  generating  waves  along  the  span 
and  essentially  pushing  entrained  water  droplets  towards  the  tunnel  walls.  After  de¬ 
position  on  the  tunnel  walls,  the  droplets  coalesce  and  drip  to  the  tunnel  hoor  forming 
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a  pool  of  suppressant.  Another  hypothesis  submits  that  droplet  transport  is  governed 
primarily  by  the  combination  of  inertial  and  fluid  dynamic  influences  instead  of  flow- 
field  dynamics  alone.  More  specifically,  the  path  of  an  entrained  particle  is  determined 
by  the  relationship  between  the  characteristic  inertial  and  flow  times  of  the  dispersed- 
phase  system.  This  relationship  ultimately  guides  the  particle  along  either  a  ballistic 
path  or  directs  the  particle  along  the  flow’s  streamlines. 

The  following  chapters  discuss  in  detail  the  theory,  methodology,  and  findings  of 
the  numerical  simulations  of  flow  past  idealized  clutter  elements  and  relate  the  results 
to  engine  nacelle  fire  suppression. 
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II.  Background  and  Theory 


As  stated  in  Section  1.1,  nacelle  fire  suppression  and  the  NGP  cover  a  wide  range 
of  topics  from  suppressant  chemistry  and  toxicity  to  delivery  and  operational 
concerns.  This  research,  however,  narrows  the  scope  of  interest  primarily  to  nacelle 
flow-held  dynamics  and  suppressant  delivery.  In  this  chapter,  discussion  is  devoted 
to  how-held  dynamics  of  a  cylinder  in  cross-how,  turbulent  how  characteristics,  and 
turbulence  modeling  techniques  such  as  Detached-Eddy  Simulation  (DES). 

2.1  Circular  Cylinder  in  Cross-Flow 

As  early  as  Leonardo  da  Vinci,  who  sketched  the  wakes  of  bluff  bodies  in  the  fifth 
century,  scientists  have  studied  the  circular  cylinder  in  cross-how.  The  phenomenon 
of  vortex  shedding  from  a  bluff  body  is  widely  known  and  has  born  significant  scrutiny 
over  a  wide  range  of  free-stream  conditions.  [5,11,18,19,28]  The  potential  solution  for 
how  past  a  circular  cylinder  predicts  symmetrical  streamlines  about  the  cylinder  and 
computes  the  how  at  the  surface  of  the  cylinder  in  polar  coordinates  as  ur  =  0  and  ue  = 
— 2U00sin  6.  In  reality,  the  potential  solution  violates  the  no-slip  condition  at  the  wall 
and  is  completely  non-physical.  Visualization  of  how  past  a  cylinder  reveals  distinct 
alternating  vortices  shed  from  the  top  and  bottom  of  the  cylinder  called  Karman 
vortex  streets.  In  1911,  von  Karman  showed  that  the  vortex  street  (Figure  2.1),  clearly 
evident  for  Reynolds  numbers  of  40  <  Rep  <  150,  is  a  stable  energy  configuration  of 
the  how.  [18]  As  the  Reynolds  number  increases  (Ren  >  105),  vortex  shedding  remains 
a  defining  feature,  but  the  vortex  streets  decay  quickly  downstream  of  the  cylinder  due 
to  turbulent  dissipation.  In  1878,  Strouhal  quantified  the  periodic  shedding  frequency 
associated  with  bluff  bodies  and  related  it  to  free-stream  velocity  and  body  diameter. 


Figure  2.1:  Depiction  of  Karman  vortex  streets,  Re  ps  160 
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The  Strouhal  number,  a  non-dimensional  parameter,  is  given  by 


St=J—  (2.1) 

where  /  is  the  frequency  of  the  vortex  shedding,  D  is  the  cylinder  diameter,  and  U 
is  the  free-stream  velocity.  The  accepted  Strouhal  number  for  Reynolds  numbers  of 
100  <  Re  <  105,  correlated  by  Rayleigh  in  1879,  is  approximately  0.2.  [18,28] 

2.2  Turbulence 

Turbulence,  often  characterized  by  high  Reynolds  numbers,  is  a  flow  regime  dis¬ 
tinguished  by  the  generation,  transport,  and  destruction  of  three-dimensional  vortex 
structures.  Smoke  issuing  from  a  smokestack  is  a  common  example  of  the  production, 
convection,  and  dissipation  of  turbulence.  [24], 

As  turbulent  eddies  move  through  the  flow,  large  amounts  of  mass,  momentum, 
and  energy  are  redistributed.  Unlike  laminar  flow  where  the  method  of  momentum 
transfer  is  molecular  diffusion,  the  dominant  mode  of  momentum  transfer  is  the  large, 
mixing  vortices.  From  an  engineering  standpoint,  turbulence  serves  to  improve  the 
mixing  of  fuel  and  oxidizer  in  combustion  reactions  and  plays  a  role  in  controlling  flow 
separation.  Energy  and  momentum  transfered  from  the  mean  flow  into  a  turbulent 
boundary  layer  serves  to  excite  the  flow  and  stave  off  separation  when  compared 
to  laminar  flow.  The  dimples  on  a  golf  ball  force  the  boundary  layer  to  become 
turbulent  in  order  to  take  advantage  of  increased  turbulent  mixing  and  prolonged 
separation.  [3,24] 

Another  defining  characteristic  of  turbulence  is  fluid  motion  over  an  expansive 
range  of  length  and  time  scales.  An  observer  watching  smoke  emanating  from  a 
stack  would  note  that  the  largest  eddies  are  on  the  order  of  the  diameter  of  the  stack. 
Likewise  in  a  turbulent  boundary  layer,  the  scale  of  the  largest  turbulent  structures  are 
on  the  order  of  the  boundary  layer  thickness,  6,  and  are  associated  with  the  integral 
length  scale.  Eddies  at  the  integral  scale  are  very  energetic  and  usually  the  direct 
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Figure  2.2:  Turbulent  boundary  layer  visualization  at  Re  =  3500  illustrating  the 

range  of  length  scales  inherent  in  turbulent  flow  [8] 

result  of  a  flow’s  turbulence  generation  mechanism.  [3]  With  respect  to  a  cylinder  in 
cross-flow,  the  integral  scale  is  on  the  order  of  the  diameter  of  the  cylinder  and  is 
observed  directly  in  a  Karman  vortex. 

Returning  to  the  smokestack  example,  an  observer  would  notice  over  time  large 
vortex  structures  breaking  down  and  giving  way  to  smaller  vortices.  Each  time  smaller 
vortices  are  born  from  a  larger  vortex,  kinetic  energy  is  transfered  to  a  smaller  scale 
in  the  flow.  The  cascade  of  energy  across  the  spectrum  of  length  and  time  scales 
is  a  fundamental  aspect  of  turbulence.  The  energy  continuously  cascades  from  the 
integral  scale  down  to  the  smallest  scale,  the  Kolmogorov  scale,  rj,  where  turbulent 
kinetic  energy  is  dissipated  directly  to  heat  via  molecular  viscosity.  [24]  Accurately 
accounting  for  the  energy  cascade  is  crucial  to  modeling  turbulence  and  computing 
its  effects.  Relative  to  the  large,  integral  scale  eddies  of  the  flow,  the  smallest  eddies 
occur  on  a  very  short  time  scale.  Kolmogorov  postulated  that  the  rate  at  which  the 
smallest  eddies  received  energy  from  the  largest  eddies  was  approximately  equal  to 
the  rate  of  dissipation  to  heat.  This  relationship  is  known  as  Kolmogorov’s  universal 
equilibrium  theory.  [29]  Computation  of  the  Kolmogorov  length  scale  based  on  the 
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universal  equilibrium  theory  supports  the  concept  that  turbulence  is  a  continuum 
phenomenon  since  rj  is  much  greater  than  the  mean  free  path  of  a  molecule,  A. 

In  summary,  turbulent  flow  is  a  seemingly  random  and  inherently  three-dimensional 
fluid  motion  characterized  by  unsteadiness  and  mixing  eddies  over  a  vast  range  of 
length  and  time  scales.  The  eddies  and  fluctuations  associated  with  the  dynamics  of 
turbulence  in  contrast  with  laminar  flow  force  definite  interactions  between  the  layers 
of  the  fluid  which  leads  to  increased  momentum  and  energy  transfer. 

2.3  Turbulence  Modeling 

Computational  Fluid  Dynamics  (CFD),  in  general,  refers  to  numerically  solv¬ 
ing  the  discretized  Navier-Stokes  equations  which  govern  the  motion  of  fluid  in  a 
continuum.  When  the  flow  physics  require  accounting  for  turbulence  in  a  numerical 
solution,  additional  computational  expenses  are  introduced.  Although  the  motion 
of  turbulence  may  seem  random,  the  fluid  motion  is  still  governed  by  the  Navier- 
Stokes  equations.  Thus,  one  method  of  computing  a  turbulent  flow  solution  is  Direct 
Numerical  Simulation  (DNS),  which  solves  the  Navier-Stokes  equations  directly  and 
resolves  the  entire  spectrum  of  turbulent  length  and  time  scales,  thus  capturing  the 
energy  cascade  discussed  in  Section  2.2.  However,  the  drawback  to  such  an  approach 
is  that  the  required  number  of  grid  points  scales  with  Re 9/4  which  leads  to  a  CPU 
time  proportional  to  Re 3  and  restricts  DNS  to  problems  with  Reynolds  numbers  on 
the  order  of  104  —  105  based  on  the  largest  current  supercomputing  capabilities.  The 
complexity  of  turbulence  usually  forces  the  adoption  of  turbulence  models  to  approx¬ 
imate  the  elements  of  turbulence  to  solve  engineering  problems  in  a  timely,  affordable 
manner.  [4] 

For  many  engineering  applications,  mean  flow  dynamics  are  the  primary  focus 
not  the  individual  fluctuations  caused  by  the  turbulent  eddies.  When  this  is  the  case, 
the  Reynolds  Averaged  Navier  Stokes  (RANS)  equations  can  be  used  to  compute  mean 
flow  quantities  by  time  averaging  the  Navier-Stokes  equations.  This  is  accomplished 
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by  substituting  an  average  and  fluctuating  component  into  the  equations  for  each  of 
the  variables.  For  instance,  velocity  becomes 


Ui  =ui+u'i 


(2.2) 


After  substituting  the  mean  plus  fluctuation  terms,  the  equations  are  averaged. 
Either  time  or  spatial  averaging  may  be  applied.  However,  time  averaging  is  appropri¬ 
ate  to  statistically  stationary  turbulence  while  spatial  averaging  is  more  appropriate 
for  homogeneous  turbulence.  Since  the  average  of  a  fluctuation  is  zero,  all  averaged 
terms  containing  a  fluctuation  are  eliminated  except  for  the  terms  containing  the  av¬ 
erage  of  the  product  of  two  fluctuations;  these  terms  are  most  certainly  non-zero.  The 
resulting  RANS  equations  are  as  follows  [4] 


»  r  dp  d  ,  . 

Mass  :  -f  —  Ips,  +  fmS)  =  0 


(2.3) 


where 


Momentum  : 


dpUi 

dt  + 


c) 

dxj 


{pUiUj 


dP  d 
dxi  dxj 


( Tij  -  pu'iU’j) 


(2.4) 


-  -  9  Q  2  9Uk A 

Tij  ^P^ij  3  ^  dXk^ 

c  _  1  (duj  duj  \ 

lj  2  \dxj  dxi ) 


(2.5) 


The  energy  equation  is  similarly  derived  using  fluctuations  and  time  averaging.  While 
R.ANS  methods  exclude  the  need  to  calculate  the  fluctuations  of  each  flow  quantity, 
R.ANS  does  however  introduce  a  new  term:  The  Reynolds  Stress  Tensor,  —Jmpuj.  The 
Reynolds  Stress  is  related  to  the  momentum  transfer  that  arises  from  the  turbulent 
fluctuations  and  introduces  a  closure  problem  to  the  R.ANS  equations.  There  are  two 
methods  for  computing  the  Reynolds  Stresses:  first  and  second  order  closures.  Second 
order  closures  are  obtained  by  deriving  exact  equations  for  the  Reynolds  stresses  by 
computing  a  second  moment.  This  introduces  additional,  higher-order,  unknown  cor- 
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relations.  Application  of  empirical  data  is  necessary  to  determine  the  new  unknowns. 
First  order  closures,  on  the  other  hand,  employ  the  eddy  viscosity  hypothesis  and 
calculate  a  turbulent  viscosity  for  the  flow.  [4] 

The  Boussinesq  eddy  viscosity  hypothesis  approximates  the  Reynolds  Stresses 
and  provides  closure  to  the  RANS  equations  by  assuming  that  the  turbulent  shear 
stress  is  related  to  the  mean  strain  rate  by  a  proportionality  factor  known  as  the  eddy 
viscosity,  pr- 


fT  = 

‘ij 


=  —pUiUj  =  2/irS, 


y 


2  duk 
xAr  w — Oij 

3  oxk 


(2.6) 


As  a  result,  the  laminar  flow  Navier-Stokes  equations  are  mathematically  iden¬ 
tical  to  RANS  formulation.  The  only  difference  is  that  the  viscosity,  p,  is  replaced 
with  the  sum  of  the  molecular  and  eddy  (turbulent)  viscosities.  In  this  research,  only 
incompressible  flows  are  of  concern,  therefore  the  terms  associated  with  the  divergence 
of  velocity,  are  neglected  in  both  Equations  (2.5)  and  (2.6). 


p  =  Pl  +  Pt  (2.7) 

The  molecular  (laminar)  viscosity  may  be  computed  using  Sutherland’s  Law 
since  it  is  a  property  of  the  fluid.  While  the  eddy  viscosity  assumption  eliminates 
the  need  to  compute  the  Reynolds  Stress  Tensor  directly,  it  does  not  eliminate  the 
requirement  to  compute  a  reasonable  estimate  of  the  turbulent  viscosity. 

Calculation  of  the  eddy  viscosity  is  accomplished  using  either  a  zero-,  one-,  or 
two-equation  model.  Zero-equation  or  algebraic  models  are  so  named  because  they  do 
not  involve  solving  a  differential  equation  model  and  instead  compute  the  eddy  viscos¬ 
ity  based  on  an  algebraic  formulation  related  to  the  velocity  profile.  A  one-equation 
model,  adds  complexity  by  solving  a  single  differential  equation  to  govern  the  trans¬ 
port,  generation,  and  destruction  of  turbulence.  The  most  elaborate  RANS  models 
employ  two  differential  equations  to  govern  turbulent  kinetic  energy  and  turbulence 
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dissipation.  Because  this  research  employs  the  Spalart-Allmaras  turbulence  model, 
attention  will  now  be  focused  on  the  development  and  theory  of  the  Spalart-Allmaras 
one-equation  turbulence  model. 

P.R.  Spalart  and  S.R.  Allmaras  developed  a  model  for  the  transport  of  a  modified 
eddy  viscosity,  z>,  based  on  Galilean  invariance,  dimensional  analysis,  and  empirical 
data.  Galilean  invariance  is  the  concept  that  fundamental  laws  of  physics  are  the 
same  for  all  reference  frames.  Spalart  and  Allmaras  developed  their  equation  for  eddy 
viscosity  term  by  term,  adding  constants  and  functions  to  calibrate  the  model,  against 
empirical  data  from  free  shear  flows,  wakes,  and  boundary  layers.  Since  Spalart  and 
Allmaras  were  developing  a  one-equation  model  focusing  on  application  to  aerody¬ 
namic  flows,  they  elected  to  abandon  the  classical  approach  of  tweaking  their  model 
for  homogeneous  turbulence  and  proceeded  directly  to  a  transport-diffusion  model 
and  calibrated  the  model  at  that  level.  The  resulting  Spalart-Allmaras  turbulence 
model  is  as  follows:  [22] 

~5t  +Uidar  =  °bl +  “  [V  '  ((^  +  7  +  C62  (Vi>)2]  ~  °wlfw  (2-8) 
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The  first  term  on  the  left  hand  side  of  Equation  (2.8)  represents  the  time  rate  of 
change  of  the  modified  turbulent  viscosity,  z>,  and  the  second  is  the  convective  term. 
Generation  of  u  is  represented  by  the  first  term  on  the  right  hand  side  and  is  followed 
by  the  turbulent  and  viscous  diffusion  terms.  The  final  term  on  the  right  hand  side 
controls  the  destruction  of  modified  turbulent  viscosity.  Also,  the  above  version  of 
the  Spalart-Allmaras  turbulence  model  does  not  include  the  trip  terms  suggested 
by  Spalart  to  allow  the  user  to  designate  transition  from  laminar  to  turbulent  flow, 
because  the  option  is  not  provided  in  Fluent™.  [9] 

To  conclude  the  discussion  of  eddy  viscosity,  one  should  note  that  although  the 
turbulent  viscosity  has  the  same  units  as  the  molecular  viscosity,  it  is  not  a  fluid 
property,  but  is  instead  representative  of  the  flow  and  is  problem  dependent. 

As  discussed  in  Section  2.2,  turbulence  is  characterized  by  a  eddies  over  a  wide 
range  of  time  and  length  scales.  The  bulk  of  the  momentum  transfer  throughout  the 
flow  is  done  by  the  largest  eddies,  which  are  more  problem  dependent.  From  the  Kol¬ 
mogorov  universal  equilibrium  theory,  the  smaller  eddies  are  generally  more  uniform 
and  easier  to  model  from  one  problem  to  the  next.  Directly  computing  the  large  ed¬ 
dies  by  using  filtered  equations  and  modeling  only  the  smallest  eddies  is  called  Large 
Eddy  Simulation  (LES).  The  governing  equations  for  Large  Eddy  Simulation  (LES) 
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are  obtained  by  filtering  the  Navier-Stokes  equations  with  a  filtering  function  based 
on  either  wave  number  or  physical  space.  [4,9]  The  resulting  filtered  equations  govern 
only  the  dynamics  of  the  eddies  larger  than  the  filter  width.  Fluent’s™  method  of 
filtering  the  governing  equations  is  implicitly  accomplished  by  using  the  cell  volume  as 
the  filtering  function.  Even  though  the  smallest  eddies  are  filtered  from  the  governing 
equations,  the  energy  cascade  to  the  smallest  eddies  cannot  be  neglected.  Subgrid 
Scale  (SGS)  models  account  for  the  energy  cascade  associated  with  eddies  smaller 
than  the  filter  width  by  modeling  the  dissipation  of  turbulent  kinetic  energy  to  heat. 
Fluent’s™  LES  SGS  models  apply  the  BoussinesqDES  is  able  to  switch  seamlessly 
between  LES  and  R.ANS  formulations  because  a  single  velocity  field  is  computed 
thereby  eliminating  discontinuities  at  the  RANS  and  LES  transition,  hypothesis  used 
in  RANS  computations  and  most  use  a  basic  mixing  length  model  to  compute  the 
eddy  viscosity.  [9]  Although  LES  allows  for  coarser  grids  and  larger  time  steps  than 
those  required  by  DNS,  the  number  of  grid  points,  N,  for  an  LES  of  a  turbulent 
boundary  layer  still  scales  with  Re2.  [14]  In  particular,  LES  requires  large  numbers  of 
cells  in  the  boundary  layer  making  it  impractical  for  many  engineering  applications. 

2.3.1  Detached-Eddy  Simulation.  In  an  effort  to  reduce  the  computational 
requirements  of  computing  turbulent  flow  while  still  maintaining  the  time-accurate 
large  eddies  that  LES  resolves,  Spalart  formulated  a  method  known  as  a  Detached- 
Eddy  Simulation  (DES).  DES  is  a  three-dimensional,  unsteady  turbulence  modeling 
technique  which  makes  use  of  both  LES  and  RANS  methods.  DES  models  switch  from 
LES  mode  to  RANS  mode  implicitly  based  on  the  grid  density  and  the  distance  to  the 
closest  wall.  When  the  model  is  operating  as  an  LES,  it  is  capable  of  resolving  eddies 
on  the  order  of  the  cell  spacing,  A.  [2,25]  The  model  carries  the  designation  Detached- 
Eddy  because  it  only  resolves  large  eddies  which  are  detached,  or  at  a  distance  from 
a  wall,  as  determined  by  the  grid  cell  size.  [25] 

In  Fluent™,  the  DES  model  uses  a  Spalart- Allmaras  turbulence  model  for  its 
RANS  mode.  While  DES  is  not  strictly  tied  to  the  Spalart-Allmaras  model,  it  does 
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offer  a  convenient  length  scale  to  use  as  the  DES  parameter.  [25]  In  order  to  trans¬ 
form  the  ordinary  R.ANS  Spalart-Allmaras  turbulence  model  into  a  DES,  a  slight 
modification  to  the  equation  is  required.  The  distance  to  the  wall,  dw,  in  the  original 
Spalart-Allmaras  transport  equation  is  replaced  by  a  quantity,  d,  which  is  the  lesser 
of  two  quantities:  the  grid  spacing,  A,  or  the  distance  to  the  wall,  dw.  [23] 

d  =  min(CDEsA,  dw)  (2.9) 

A  =  max  (Ax,  Ay,  A z) 

The  constant,  Cdes ,  is  typically  set  to  0.65  based  on  calibration  in  isotropic  turbu¬ 
lence. 

With  this  modification  to  the  Spalart-Allmaras  length  scale,  the  model  switches 
automatically  from  a  RANS  formulation  near  the  wall  to  an  LES/SGS  away  from  the 
wall.  When  using  a  pure  RANS  model  throughout  the  domain,  the  eddy  viscosity 
is  designed  to  drop  off  to  zero  in  the  free-stream  where  no  turbulence  is  modeled. 
By  altering  the  length  scale  of  the  Spalart-Allmaras  turbulence  model,  the  turbulent 
production  and  destruction  terms  in  Equation  (2.8)  cannot  become  zero  in  the  free- 
stream  thereby  allowing  the  generation,  transport,  and  dissipation  of  eddy  viscosity 
on  the  subgrid  scale.  In  the  LES  region,  turbulent  structures  smaller  than  the  grid 
cell  size  are  modeled  by  the  Spalart-Allmaras  turbulence  model  functioning  as  an 
SGS  eddy  viscosity  model  while  larger  eddies  in  the  flow  are  directly  resolved.  DES 
is  able  to  switch  seamlessly  between  LES  and  RANS  formulations  because  a  single 
velocity  field  is  computed  thereby  eliminating  discontinuities  at  the  RANS  and  LES 
transition. 

The  behavior  of  the  model  is  controlled  by  the  user  through  grid  design  by 
controlling  A.  For  example,  the  user  should  build  a  grid  where  A  6  >  dw  to  trigger 
RANS  behavior.  [25]  Conversely,  in  regions  where  an  LES  response  is  desired,  the  grid 
density  should  correspond  to  d  =  CdesA. 


15 


Grid  refinement  of  a  DES  grid  is  often  challenging  because  A,  the  local  grid 
spacing,  is  the  largest  value  of  the  spacing  in  the  x,  y,  and  z  directions.  Consequently, 
refinement  in  only  one  or  two  directions  would  be  wasted;  isotropic  grid  density  is  the 
goal  in  the  LES  region.  On  the  other  hand,  the  severe  anisotropy  (high  aspect  ratio) 
of  boundary  layer  grids  is  the  mechanism  for  triggering  a  RANS  response  near  a  wall. 

Once  more,  the  only  parameter  that  varies  to  alter  the  LES/RANS  state  is  d. 
In  engineering  situations  where  resolution  of  the  boundary  layer  eddies  is  not  criti¬ 
cal,  DES  offers  a  method  for  directly  computing  the  large,  turbulent  flow  structures 
dominating  massively  separated  flows  at  an  affordable  computational  cost. 

2.3.2  Free-stream  Turbulence  Intensity  for  Spalart-Allmaras  Based  DES. 

The  turbulence  intensity  of  a  flow,  which  is  defined  as 

/=—  (2.10) 

^ avg 

where  u'  is  the  RMS  of  the  velocity  fluctuations  and  uavg  is  the  mean  velocity,  can 
significantly  impact  separation  behavior.  However,  when  using  a  DES  model  based 
on  the  Spalart-Allmaras  turbulence  model,  specification  of  free-stream  turbulence 
intensity  has  little  impact  on  the  outcome  of  the  solution.  In  the  absence  of  large 
scale  eddies  for  direct  evaluation,  the  Spalart-Allmaras  SGS  model  estimates  the 
effect  of  turbulence  exclusively  with  z>,  the  modified  turbulent  viscosity.  As  a  result, 
setting  a  free-stream  turbulence  intensity  at  the  inlet  only  serves  to  alter  the  value  of 
the  modified  turbulent  viscosity  at  the  boundary.  Furthermore,  the  Spalart-Allmaras 
turbulence  model  is  insensitive  to  high  values  of  z>  in  the  free-stream  and  damps 
disturbances  quickly  in  the  absence  of  turbulence  generation.  In  the  current  effort, 
free-stream  turbulence  intensity  would  be  an  interesting  variable  to  examine,  but 
the  limitations  of  the  Spalart-Allmaras  SGS  model  prevent  a  meaningful  analysis  of 
free-stream  turbulence  intensity  when  using  DES.  [4,22] 
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Figure  2.3:  DES  Grid  Regions;  ER  =  Euler  Region,  RR  =  RANS  Region,  FR  = 
Focus  Region,  DR  =  Departure  Region 


2.3.3  DES  Grids.  Since  DES  combines  RANS  and  LES  methods  into  a 
single  model,  the  complications  associated  with  generating  RANS  and  LES  grids 
individually  are  also  present  in  a  DES  grid  generation.  As  discussed  in  Section  2.3.1, 
grid  density  controls  the  path  a  DES  travels:  RANS  or  LES.  Before  designing  a  DES 
grid,  a  user  should  be  intimately  familiar  with  the  expected  flow-held  in  each  region 
of  the  grid. 

A  DES  grid  is  divided  into  three  Super- Regions:  the  Euler  (ER),  RANS  (RR), 
and  LES  (LR)  regions.  The  Enler  Region  of  a  DES  mesh  is  typically  located  up- 
stream  of  any  turbulence  generation  and  has  cell  spacing  dictated  only  by  geometry 
and  shocks.  Generally,  coarse,  isotropic  cell  spacing  in  all  three  directions  prevails.  [21] 
The  R  ANS  Region  of  the  grid  contains  the  how’s  boundary  layer  and  is  itself  comprised 
of  the  RANS  Viscous  and  Outer  Regions.  The  RANS  Viscous  Region  under  DES  is 
identical  to  the  near  wall  region  of  a  purely  RANS  grid  and  is  responsible  for  capturing 
the  viscous  sublayer,  buffer  layer  and  the  log  layer.  When  using  the  Spalart-Allmaras 
turbulence  model  in  the  RANS  Region,  Spalart  recommends  the  spacing  of  the  hrst 
cell  in  the  wall-normal  direction  be  at  least  A y+  =  2  or  less.  [21]  Fluent™  suggests  tar¬ 
geting  a  A y+  of  one  for  the  Spalart-Allmaras  model.  [9]  Additionally,  the  wall-normal 
growth  rate  should  be  on  the  order  of  1.25  or  less  to  accurately  capture  the  log  layer. 
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In  the  wall-parallel  direction,  grid  spacing  is  dictated  by  standard  RANS  methods 
and  therefore  shocks  and  geometry  drive  the  requirements  for  the  Aa;+  grid  refine¬ 
ment.  [21]  In  the  RANS  Outer  Region,  standard  RANS  grid  spacing  practices  prevail. 
The  cause  for  concern  in  this  region  is  the  free-stream/boundary  layer  interface.  In 
the  free-stream,  eddy  viscosity  is  near  zero  and  may  rapidly  approach  zero  at  the 
edge  of  the  RANS  region.  For  an  accurate  model,  the  boundary  layer  solution  must 
detect  a  slope  discontinuity  in  the  eddy  viscosity  across  the  free-stream/boundary 
layer  interface.  As  a  result,  a  mesh’s  RANS  Region  is  often  extended  well  past  the 
boundary  layer  height,  5.  The  solution  is  far  more  forgiving  when  5  is  overestimated 
rather  than  underestimated.  [21] 

The  LES  Region  is  subdivided  into  two  regions:  the  Focus  (FR)  and  Departure 
(DR)  Regions.  As  the  name  implies,  the  Focus  Region  is  the  area  of  the  model  where 
accurate  resolution  of  the  detached-eddies  is  desired.  The  objective  cell  spacing,  A0, 
in  the  Focus  Region  quantifies  the  spatial  resolution  of  a  DES.  Grid  spacing  in  the  FR 
should  be  isotropic  since  A  =  max( Ax,  Ay,  A z)  defines  the  local  grid  density  and 
serves  to  filter  out  turbulent  eddies  on  the  order  of  the  local  grid  cell  size.  Once  more, 
the  user  should  remember  that  refining  the  FR  mesh  in  only  one  or  two  directions 
wastes  computational  resources  (Section  2.3.1).  Finally,  the  Departure  Region  is 
found  adjacent  to  the  FR  and  transitions  the  LES  grid  spacing  of  A0  to  coarser 
cell  spacing  often  on  the  order  of  Euler  spacing  near  grid  boundaries.  Essentially, 
the  DR  serves  as  a  good  numerical  neighbor  to  the  Focus  Region  and  only  imposes 
a  moderate  computational  cost  on  the  solver.  Results  in  this  region  of  the  mesh 
should  not  be  accepted  as  spatially  accurate  since  the  motion  of  the  resolved  eddies 
is  increasingly  dominated  by  eddy  viscosity  as  the  grid  density  coarsens.  Spalart 
suggests  the  following  rule  of  thumb  for  determining  how  large  to  make  the  Focus 
Region  and  where  to  start  the  Departure  Region.  “Can  a  particle  return  from  this 
[DR]  point  to  very  near  the  body?  Is  there  flow  reversal?”  [21]  If  the  answer  to  these 
questions  is  yes,  then  regions  containing  such  phenomena  should  be  included  in  the 
FR, 
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2.3.4  DES  Time  Step.  Selection  of  a  suitable  time  step  for  a  Detached- 
Eddy  Simulation  (DES)  is  based  on  the  the  grid  density  of  the  Focus  Region  because 
the  FR  possesses  the  finest  spatial  resolution  where  unsteady  turbulent  structures  are 
resolved.  None  of  the  other  regions  will  likely  present  frequency  content  higher  than 
the  FR.  Assuming  that  the  SGS  model  is  well  calibrated,  then  the  smallest,  potentially 
active  eddies  would  have  a  wavelength  A  ~  5A.  [21]  Although  these  eddies  are  resolved, 
they  are  not  particularly  spatially  accurate  because  their  motion  is  governed  by  eddy 
viscosity  and  the  cascade  of  energy  to  smaller  eddies  is  modeled  instead  of  directly 
computed.  Given  these  conditions,  Spalart  recommends  five  time  steps  per  period  for 
a  local  Courant  Number  (CFL)  of  1.  This  means  that  a  wave  (eddy)  travels  across 
one  cell  for  every  timestep.  Thus,  five  samples  per  period  will  give  a  rudimentary 
sampling  of  the  eddy.  Based  on  these  assumptions,  a  good  starting  point  for  selecting 
the  solver  time  step,  At  is  to  apply  the  formula 

A  t=A-  (2.11) 

U  max 

where  A0  is  the  targeted  FR  cell  spacing  and  Umax  is  the  maximum  estimated  velocity 
encountered  in  the  FR.  Spalart  reminds  the  user  that  this  formulation  is  only  a  start¬ 
ing  point,  and  space/time  balancing  is  encouraged  since  not  all  differencing  schemes 
behave  alike.  [21] 

2.3.5  Grid  Convergence.  Ideally  for  a  DES  solution,  simulations  are  con¬ 
ducted  at  cell  spacings  of  A0  and  A0/2.  Of  course,  for  some  grids,  halving  the  grid 
spacing  in  the  FR  can  present  an  impossible  computational  burden  so  at  the  very 
least,  the  number  of  grid  cells  between  the  nominal  grid  and  the  convergence  test  grid 
should  increase  by  a  factor  of  the  y/2.  [21]  Due  to  the  implicit  relationship  between 
space  and  time  for  a  LES,  both  the  grid  spacing  and  the  time  steps  size  must  be 
refined  in  order  to  conduct  a  complete  grid  convergence  study. 
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2-4  Droplet  Transport 


2-4-1  Inertial  Analysis.  A  simple,  first-order  estimate  of  dispersed-phase 
particle  motion  assumes  that  a  droplet’s  path  is  governed  exclusively  by  inertia  and 
the  local  velocity  vector  of  a  carrier  fluid.  For  this  analysis,  the  vertical  component 
of  droplet  drag  is  negligible,  and  an  entrained  droplet  is  assumed  to  move  at  a  speed 
equal  to  the  magnitude  of  the  local  velocity  of  the  carrier  fluid.  However,  the  entrained 
droplet  does  not  necessarily  track  the  carrier  fluid’s  streamlines  exactly.  Instead,  the 
droplet  traverses  a  path  described  by  the  vector  sum  of  the  carrier  velocity  and  the 
velocity  imposed  by  gravity.  When  this  is  the  case,  a  droplet  entrained  in  a  fast 
moving  carrier  fluid  could  pass  through  the  clutter  array  before  gravity  causes  the 
particle  to  fall  into  one  of  the  cylinders.  Conversely,  in  slower  flow-fields,  the  time 
required  to  pass  the  array  may  be  significantly  longer  or  on  the  order  of  the  time 
required  to  fall  into  one  of  the  clutter  elements. 


2-4-2  Stokes  Number.  Further  analysis  of  droplet  velocity  and  inertia  in¬ 
volves  computing  the  Stokes  number  throughout  the  flow-held.  The  Stokes  number 
is  defined  as  the  ratio  of  the  response  time  of  a  dispersed  particle  to  the  response 
time  of  the  system.  Essentially,  the  Stokes  number  is  a  non-dimensional  parameter 
for  characterizing  the  ability  of  entrained  particles  to  negotiate  obstacles.  The  Stokes 
number  is  given  by 

Stk  =  —  (2.12) 

ts 

where  the  characteristic  time  of  the  droplet,  r^,  is  given  by 


Td  = 


(2.13) 


and  pd  is  the  density  of  the  droplet,  dd  is  the  diameter  of  the  droplet,  and  pc  is  the 
viscosity  of  the  carrier  fluid.  Finally  ts,  the  system  response  time,  is  computed  based 
on  a  characteristic  length,  LSl  and  a  characteristic  velocity,  Vs.  With  respect  to  the 
cylinder  array,  the  characteristic  length  is  the  cylinder  diameter  because  it  poses  the 
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most  conservative  estimate  for  the  response  time  of  the  system.  For  Stokes  numbers 
Stk  <  1.0,  the  dispersed  phase  particles  follow  the  streamlines  of  the  carrier  fluid 
closely.  On  the  other  hand,  for  Stk  >1.0  the  path  of  the  particles  is  increasingly 
governed  by  inertia  and  not  the  carrier  fluid’s  streamlines.  The  smaller  the  Stokes 
number,  the  better  a  droplet  will  be  able  to  negotiate  an  obstacle.  [9] 

2.5  CFD  Cylinder  Research 

A  single  circular  cylinder  in  cross-flow  has  been  studied  by  a  multitude  of  CFD 
RANS  models.  [2,25,26,28,29]  However,  with  the  advent  of  DES,  many  investigators 
have  chosen  to  exercise  DES  on  the  circular  cylinder  problem  and  validate  it  as  a  viable 
option  for  computing  bluff  body  flows.  Travin  et  al.  modeled  laminar  and  turbulent 
separated  flow  past  a  circular  cylinder  using  DES  and  found  good  agreement  between 
the  model  and  experimental  data  for  drag,  shedding  frequency,  pressure,  and  skin 
friction.  [25]  In  another  study,  Vatsa  et  al.  compared  low  speed  flow  past  a  circular 
cylinder  with  Reynolds  number  ranging  from  5xl04  to  1.4xl05  using  NASA  Langley’s 
production  TLNS3D  code  for  2-D  and  3-D  URANS  and  DES  models.  Vatsa’s  group 
found  that  3-D  DES  matched  well  with  the  experimental  data  and  work  of  Travin 
et  al.  and  judged  3-D  DES  an  acceptable  method  for  predicting  highly  separated 
flows.  [26] 

Unlike  the  single  circular  cylinder  in  cross-flow,  an  array  of  cylinders  has  borne 
little  scrutiny  by  experimentalists  and  CFD  practitioners  alike.  Most  existing  research 
in  the  realm  of  flow  past  an  array  of  cylinders  was  accomplished  for  the  purpose 
of  determining  and  predicting  heat  transfer  and  vibration  in  heat  exchanger  tube 
bundles.  Using  a  2-D  RANS  method,  Watterson  et  al.  modeled  the  flow  around 
a  staggered  array  of  cylinders  in  turbulent,  low-speed  cross  flow  (Ren  =  21,000). 
The  numerical  scheme  of  Watterson  et  al.  [27]  solved  the  governing  Navier-Stokes 
equations  via  a  pressure  correction  method  and  applied  the  low-Reynolds  number 
k-e  turbulence  model  of  Yang  and  Shih.  Wind  tunnel  results  were  compared  to  the 
numerical  predictions  of  Watterson’s  group  and  achieved  only  marginal  agreement 
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with  experimental  trends.  [27]  Seeking  a  higher  fidelity  model,  Hassan  and  Barsamian 
explored  the  staggered  tube  bundle  problem  but  instead  employed  a  LES  technique 
allowing  spatial  and  temporal  resolution  of  the  large-scale  turbulent  eddies  associated 
with  the  tube  wakes  at  Re.o  —  21,  700.  The  LES  results  correlated  very  well  with  the 
same  experimental  data  used  by  Watterson  et  al.,  and  captured  the  high  transverse 
turbulence  intensities  and  the  generation  of  span-wise  vorticity  aspects  invisible  to  a 
2-D  RANS  solution.  [12] 
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III.  Methodology 


0  taming  a  numerical  solution  for  characterizing  the  flow-held  dynamics  inside  an 
array  of  cylinders  was  an  iterative  and  nearly  simultaneous  undertaking.  The 
methodology  of  grid  generation,  numerical  simulation,  and  data  reduction  were  far 
from  serial  processes  because  each  aspect  was  implicitly  driven  by  the  others.  This 
chapter  discusses  each  process  in  the  solution  and  considers  the  links  between  each 
step. 

3. 1  Clutter  Geometry 

The  clutter  geometry  for  this  research  is  patterned  after  the  array  of  circular 
cylinders  used  by  Disimile  et  al.  [7]  The  principal  difference  between  the  two  clutter 
packages  is  scale.  The  full-size  wind  tunnel  setup  used  by  Disimile  incorporated 
cylinders  with  diameters  D  =  50.8  mm  (2  in.)  in  a  tunnel  with  a  610  x  914  mm 
(24  x  36  in.)  cross-section.  The  clutter  package  for  this  research,  however,  was  scaled 
down  by  a  factor  of  two  producing  cylinder  diameters  of  D  —  25.4  mm  (1  in.)  and 
a  tunnel  cross-section  of  304.8  x  304.8  mm  (12  x  12  in.).  This  was  done  in  order 
to  match  the  experimental  setup  of  an  ongoing  parallel  research  effort.  Figure  3.1 
represents  the  clutter  geometry  used  in  this  research. 

3.2  Test  Matrix 

Due  to  the  significant  computational  time  required  to  perform  Detached-Eddy 
Simulation  (DES),  the  scope  of  the  investigation  was  narrowed  to  ensure  data  collec¬ 
tion  was  completed  in  a  timely  manner.  The  work  of  Disimile  et  al.  [7]  indicated  that 
low-speed  (0.5—2  m/s)  flow  coupled  with  tight  stream-wise  cylinder  spacing  (0.25D) 
posed  a  significant  obstacle  to  water  droplet  transport  past  an  array  of  circular  cylin¬ 
ders  acting  as  nacelle  clutter.  For  this  research,  the  cylinder  spacing  was  fixed  at 
2.0  D. 

The  flow  dynamics  for  an  array  of  circular  cylinders  is  not  correlated  to  Reynolds 
number  like  it  is  for  a  single  cylinder.  However,  matching  the  Reynolds  number 
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Figure  3.1:  Idealized  Engine  Nacelle  Clutter  for  CFD  Model 

of  Disimile’s  experiment  to  the  CFD  model  was  considered  good  practice 
context,  the  Reynolds  number  of  the  cylinder  array  was  defined  as 

V 

where  U  is  the  free-stream  velocity,  D  is  the  diameter  of  a  cylinder,  and  v  is  the 
kinematic  viscosity.  Because  this  research  seeks  to  characterize  mechanisms  hindering 
the  dispersement  and  transport  of  fire  suppressant  droplets,  the  low-speed  realm  of 
Disimile’s  research  was  targeted.  To  maintain  Reynolds  number  consistency  between 
the  experiment  of  Disimile  et  al.  and  the  numerical  model,  a  velocity  of  1  m/s  was 
used  for  the  lower  bound  to  match  the  0.5  m/s  of  Disimile  et  al.  The  5  m/s  case, 
corresponding  to  Disimile’s  2.5  m/s  case  based  on  Reynolds  number  scaling  provides 
another  test  point  to  explore  and  replicate  the  trends  observed  by  Disimile  et  al. 
without  reaching  the  upper  bound  of  suppressant  transport  noted  by  Disimile. 

Table  3.1  compiles  the  six  test  cases  examined  in  this  research  used  to  search 
for  mechanisms  of  droplet  transport  inside  the  clutter  array.  Instead  of  examining 
the  test  matrix  one  case  at  a  time,  the  investigation  progressed  as  a  simultaneous 


In  this 

(3.1) 
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Table  3.1:  Test  Matrix  of  Numerical  Simulations 


Single  Cylinder 

Cylinder  Array  (9.9M  Cells) 

Cylinder  Array  (14M  Cells) 

A0  =  0.2 

A0  =  0.2 

A0  =  0.1 

Uoo  =  1  m/s 

Uoo  =  1  m/s 

Uoo  —  1  rn/s 

At  =  0.01  s 

At  =  0.01  s 

At  =  0.01  s 

Uoo  =  5  m/s 

Uoo  =  5  m/s 

Uoo  —  5  m/s 

At  =  0.005  s 

At  =  0.005  s 

At  =  0.005  s 

and  iterative  process.  To  start,  2-D  single  cylinder  and  the  cylinder  array  grids  were 
generated  for  computing  RANS  solutions  to  test  initial  estimates  of  boundary  layer 
near-wall  spacing.  Between  each  simulation,  the  grids  were  tweaked  until  the  wall  y+ 
proved  satisfactory  and  the  boundary  layer  profile  was  adequately  captured.  Next, 
the  2-D  grids  were  extended  into  three  dimensions,  and  RANS  computations  were 
again  performed  to  verify  wall  y+  and  the  boundary  layer  profile  and  height.  With 
satisfactory  RANS  results  in  hand,  Detached-Eddy  Simulations  were  performed  on 
the  single  cylinder  grid  to  validate  the  DES  model.  Upon  achieving  satisfactory  DES 
results  for  a  single  cylinder,  Detached-Eddy  Simulations  were  accomplished  on  the 
cylinder  array  grids.  The  following  sections  provide  deeper  insight  into  grid  generation 
and  the  solution  process.  The  grid  spacing  and  dimensions  discussed  are  the  final 
values  after  iterating  through  several  versions  of  each  grid. 

3.3  RANS  Grid  Generation 

Grid  generation  was  accomplished  using  the  commercial  software,  Gridgen™  to 
generate  2-D  and  3-D  unstructured  grids  for  DES.  Developing  node  connectivity  for 
an  array  of  cylinders  with  2.0 D  stream-wise  spacing  was  the  first  task.  After  lay¬ 
ing  out  the  cylinder  and  tunnel  wall  locations,  each  cylinder  was  wrapped  with  a 
structured  block  containing  hexahedral  (quadrilateral  in  2-D)  cells  which  employed 
RANS  spacing  methods.  Planning  for  the  Spalart-Allmaras  turbulence  model,  initial 
wall-normal  spacing  was  set  at  7.25xl0-4  for  the  hyperbolically-generated,  structured 
boundary  layer  which  grew  at  a  rate  of  1.125  for  fifteen  cells  and  then  switched  to 
a  growth  rate  of  1.5  until  the  cells  reached  approximately  unit  aspect  ratio.  For 
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an  exclusively  RANS  computation,  the  span-wise  spacing  was  not  critical  and  was 
therefore  set  equal  to  the  nominal  size  of  the  tetrahedral  (triangular  in  2-D)  cells 
comprising  the  majority  of  the  domain.  Generation  of  the  unstructured,  tetrahedral 
(triangular  in  2-D)  cells  was  accomplished  using  Delaunay  triangulation  in  Gridgen™. 
In  order  to  save  on  computational  time,  the  boundary  layers  on  the  tunnel  walls  were 
not  modeled.  The  coarser  wall-normal  grid  density  associated  with  inviscid  walls  still 
captured  pressure  gradients  and  reflections  from  the  wall,  but  significantly  reduced 
data  storage  and  computational  requirements. 


3.4  DES  Grid  Generation 

The  RANS  grid  described  in  Section  3.3  eventually  evolved  into  a  DES  grid  after 
considering  factors  pertinent  only  to  DES.  Since  the  user  indirectly  controls  when  the 
solver  switches  from  RANS  to  LES,  properly  estimating  the  boundary  layer  is  crucial 
to  determining  appropriate  grid  spacing.  Prandtl’s  formulation  for  the  thickness  of  a 
turbulent  boundary  layer  on  a  flat  plate  was  applied. 


£99  — 


0.37a; 


(3.2) 


where  a;  is  a  characteristic  length  and  Rex  is  the  Reynolds  number  based  on  x.  This 
estimation  gives  hgg  =  0.06  in.  The  quarter-circumference  of  the  cylinder  was  desig¬ 
nated  as  the  characteristic  length  because  the  critical  boundary  layer  height  occurs 
at  a  location  near  the  top  and  bottom  of  the  cylinder.  Notably,  Prandtl’s  estimation 
assumes  a  flow  with  a  zero  pressure  gradient,  which  is  of  course  not  the  case  for  a 
cylinder  in  cross-flow.  However,  the  estimate  does  gives  a  reasonable  approximation 
for  the  upper  bound  of  the  RANS  region  when  designing  the  grid. 

Due  to  the  relatively  small  diameter  of  the  cylinders,  the  cell  spacing  circum¬ 
ferentially  along  the  cylinder  wall  competed  with  the  wall  normal  spacing  for  the 
boundary  layer  mesh.  In  order  to  adequately  capture  flow  gradients  resulting  from 
the  curvature  of  the  cylinders,  an  increased  number  of  grid  points  were  placed  along 
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Region  (ER) 


Focus  Region  (FR)  Departure 

Region  (DR) 


Figure  3.2:  DES  Cylinder  Array,  9.9  million  cells 

the  wall  which  lowered  the  aspect  ratio  of  the  cells  in  the  RANS  region.  The  lower 
aspect  ratio  reduced  the  anisotropy  of  the  RANS  cells  which  in  turn  brought  the  LES 
transition  closer  to  the  cylinder  surface.  Another  competing  factor  was  the  hexahe- 
dral  to  tetrahedron  transition.  Tetrahedrons  were  intended  for  the  domain  outside  of 
the  boundary  layer.  To  reduce  high  skewness  and  the  associated  grid  metric  errors, 
hexahedrons  with  aspect  ratios  of  order  one  were  desired  at  the  boundary  between 
hexahedral  and  tetrahedral  cells.  Resolution  of  these  competing  factors  resulted  in 
boundary  layer  hexahedrons  growing  significantly  past  £99  and  into  the  LES  region, 
which  forced  some  degree  of  anisotropy  to  be  accepted  for  the  LES. 

The  Focus  Region  (FR),  defined  in  Section  2.3.3  and  delineated  in  Figure  3.2, 
is  characterized  by  fine,  isotropic  cell  spacing  designed  for  the  LES  branch  of  DES. 
On  the  cylinder  array  mesh,  the  FR  extended  two  diameters  (2 D)  upstream  of  the 
first  row  of  cylinders  and  six  diameters  (6 D)  downstream  of  the  last  row  of  cylinders. 
The  same  FR  proportions  were  also  applied  to  the  single  cylinder  grid.  Selection  of 
the  FR  objective  grid  density,  A0,  is  discussed  in  Section  3.4.1. 


3-4-1  Space/Time  Coupling.  Arriving  at  a  A0  suitable  for  capturing  the 
necessary  flow  physics  without  exceeding  the  capabilities  of  the  available  computing 
resources  required  balancing  spatial  and  temporal  concerns.  Based  on  an  expected 
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Figure  3.3:  RANS-LES  Transition,  ynorm  =  0.065  for  the  14  million  cell  grid 

Strouhal  number  of  0.2,  the  frequency  of  the  bulk  shedding  is  on  the  order  of  8  Hz 
for  the  1  m/s  case  and  on  the  order  of  40  Hz  for  the  5  m/s  case.  Spalart’s  guidance 
for  estimating  the  time  step  given  by  Equation  (2.11)  suggests  ambitious  timesteps  of 
At  =  0.04  s  for  the  1  m/s  inlet  case  and  At  =  0.015  s  for  the  5  m/s  inlet  case  based 
on  A0  =  0.2.  Determining  the  physical  timestep  of  a  CFD  solver  essentially  sets  the 
sampling  rate  thereby  presetting  the  highest  discernible  spectral  content  of  the  data. 
The  Nyquist  sampling  rate,  given  by  2/  where  /  is  the  frequency  of  a  desired  signal, 
dictates  the  minimum  sampling  rate  required  to  resolve  /.  Consequently,  utilization 
of  more  conservative  timesteps,  At  =  0.01  and  At  =  0.005,  provided  better  resolution 
of  the  expected  flow  physics.  Timesteps  such  as  these  yield  sampling  rates  equivalent 
to  100  Hz  for  the  1  m/s  case  and  200  Hz  for  the  5  m/s  case. 

In  the  end,  the  grid  density  for  the  single  cylinder  and  the  9.9  million  cell 
grids  was  A0  =  0.2,  a  compromise  between  sampling  the  flow  physics  and  conserving 
computational  time.  The  fine,  14  million  cell  grid,  which  halved  Ao  to  0.1  in  the 
Focus  Region,  was  constructed  for  testing  grid  convergence  to  ensure  the  numerical 
solution  was  independent  of  grid  density.  [21] 
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Figure  3.3  illustrates  the  RANS-DES  transition  at  yn0rm  —  0.065,  based  on  the 
14  million  cell  grid  where  Ao  =  0.1.  Since  Ay  (wall-normal)  in  the  boundary  layer  is 
very  small  compared  to  Ax  (circumferential)  and  A^  (span- wise),  the  circumferential 
and  span-wise  spacing  drive  the  LES  transition  point.  In  the  case  of  the  14  million 
cell  grid,  Ax  =  0.0877  and  Az  =  0.1  (based  on  A0),  and  the  local  grid  density  is 
computed  by  A  =  max( Ax,  Ay,  A z)  =  0.1.  Thus,  the  RANS  to  LES  transition 
point  occurs  when  Cdes A  =  0.065  <  dw. 

3.5  Fluent  Setup 

This  section  discusses  the  basic  setup  for  the  CFD  solver,  Fluent™  6.2.16 ,  for 
both  RANS  and  DES. 

Once  the  grid  design  was  completed  in  Gridgen™,  it  was  exported  to  a  Fluent™  case 
hie.  Case  hies  contain  the  grid  geometry,  grid  points,  boundary  conditions,  and  solver 
information.  After  starting  the  appropriate  version  of  Fluent™,  serial  2-D/3-D  dou¬ 
ble  precision,  a  grid  check  was  conducted  followed  by  a  series  of  smooth/swap  cycles 
to  reduce  skewness  and  improve  grid  quality.  To  ensure  the  solver  interpreted  the 
correct  dimensions  of  the  grid  geometry,  the  grid  was  scaled  by  inputting  the  units 
used  to  create  the  grid.  Following  these  grid  operations,  the  case  hie  was  saved  and 
Fluent™  was  restarted  using  the  parallel  mode  to  enable  grid  partitioning. 

Depending  on  the  availability  of  processors  on  the  Beowulf  cluster,  the  nominal 
9.9  million  cell  grid  was  partitioned  into  twelve  to  sixteen  partitions  using  the  Metis 
method.  The  14  million  cell  convergence  grid  was  partitioned  for  twenty  processors 
and  the  single  cylinder  case  with  2.2  million  cells  was  partitioned  for  four  processors. 

Even  though  the  Navier-Stokes  equations  comprise  a  set  of  nonlinear,  coupled 
equations,  the  Segregated  Solver  in  Fluent™  solves  the  governing  equations  sequen¬ 
tially  and  segregated  from  one  another.  In  the  incompressible,  ideal  gas  realm,  the 
absence  of  strong  gradients  resulting  from  shocks  allows  faster  convergence  with  an 
iterative  sequential  solver  than  with  a  simultaneous,  coupled  solver.  Fluent’s™  pro- 
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cedure  for  sequentially  solving  the  governing  equations  is  as  follows.  First,  the  fluid 
properties  are  updated  throughout  the  domain;  if  the  solver  is  on  the  first  iteration, 
the  initialized  fluid  properties  are  used.  Next,  the  velocity  field  is  updated  by  solving 
the  momentum  equations  for  u,  v,  and  w  using  the  current  values  of  pressure  and 
mass  flux  at  cell  faces.  Because  the  fluid  is  set  to  ideal,  incompressible  gas,  den¬ 
sity  is  not  explicitly  solved  as  a  conserved  variable.  To  ensure  continuity  is  satisfied, 
Fluent™  uses  a  pressure-correction  method  to  adjusted  the  velocity  held  to  balance 
mass-fluxes  at  cell  faces.  With  the  corrected  velocity  held,  scalar  equations  such  as 
energy  and  turbulence  are  solved.  Finally,  a  convergence  check  is  accomplished;  if  the 
solution  is  not  converged,  the  solver  continues  to  iterate.  [9] 

When  employing  an  ideal  gas  assumption  with  M  <1,  the  changes  in  pressure 
throughout  the  solution  are  very  small.  This  increases  the  solution’s  susceptibility  to 
round-off  error.  To  combat  the  round-off  error,  the  Operating  Pressure  was  set  to  the 
free-stream  pressure  of  101,325  Pa,  allowing  Fluent™  to  solve  for  gage  pressure,  the 
difference  between  local  static  pressure  and  operating  (reference)  static  pressure.  [9] 

Often  overlooked,  solver  boundary  conditions  deserve  consideration  in  this  sec¬ 
tion.  The  computational  domain  was  essentially  a  wind-tunnel  model:  a  rectangular 
prism  with  inlet  and  outlet  planes  at  opposite  ends  of  the  tube.  Fluent’s™  Veloc¬ 
ity  Inlet  boundary  condition,  specifically  designed  for  low-speed  incompressible  hows, 
was  invoked  at  the  inlet  plane.  At  the  Velocity  Inlet  boundary  condition,  scalar 
quantities  such  as  velocity  components,  temperature,  and  modified  turbulent  viscos¬ 
ity  were  specihed  while  stagnation  values  were  allowed  to  hoat  and  normalize  to  the 
velocity  distribution.  For  a  cell  adjacent  to  a  Velocity  Inlet,  Fluent™  computed  the 
momentum  and  energy  fluxes  across  the  cell  face  normal  to  the  boundary  with  the 
user’s  velocity  inputs.  [9]  On  the  outlet  plane,  an  Outflow  boundary  condition  was 
applied.  Outflow  boundaries  are  suitable  when  the  velocity  and  pressure  at  the  outlet 
are  unknown  and  downstream  conditions  do  not  significantly  impact  the  upstream 
how  dynamics.  Fluent™  handled  the  Outflow  boundaries  with  a  simple  extrapola¬ 
tion  from  the  grid  interior.  Finally,  the  walls  of  the  tunnel  were  modeled  as  slip  walls. 
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This  was  accomplished  in  Fluent  by  specifying  zero  shear  stress  under  the  Momentum 
panel  for  a  Wall  boundary  condition.  The  Wall  boundary  condition  at  the  cylinder 
surface  used  the  default  No-Slip  condition  for  the  wall  shear  thereby  setting  the  mod¬ 
ified  turbulent  viscosity,  z>,  equal  to  zero.  Note  that  since  the  grid  is  sufficiently  fine 
to  resolve  the  laminar  sublayer  ( y+  =  1  for  Spalart-Allmaras  turbulence  model),  the 
shear  stress  computation  at  the  wall  used  the  inner  law  relationship 

u  yv* 
v*  v 

where  v*  is  the  friction  velocity  given  by 

(3.4) 

and  u  is  the  mean  velocity  parallel  to  the  wall,  y  is  the  normal  distance  from  the 
wall,  v  is  the  kinematic  viscosity,  and  tw  is  the  shear  stress  at  the  wall.  The  thermal 
condition  for  both  No-Slip  and  Slip  Walls  was  set  to  zero  heat  flux  at  the  cylinder 
wall.  [9,28] 

3.5.1  RANS  Specific  Setup.  For  RANS  computations,  the  turbulence  model 
was  set  to  Spalart-Allmaras  on  the  Viscous  panel  and  all  default  values  were  accepted. 
The  Segregated  Solver  was  set  to  use  a  second-order  upwind  discretization  for  mo¬ 
mentum,  energy,  and  modified  turbulence  viscosity  achieved  by  applying  a  Taylor 
series  expansion  of  these  cell-centered  quantities  about  the  cell  centroid.  Finally,  exe¬ 
cution  of  a  R.ANS  computation  was  accomplished  by  initializing  the  entire  domain  to 
the  inlet  (free-stream)  conditions  and  iterating  until  convergence  was  achieved.  The 
convergence  criteria  for  continuity,  velocity,  and  turbulent  viscosity  was  set  at  three 
orders  of  magnitude.  However,  lift  and  drag  on  the  cylinder  array  were  often  used  to 
quantify  convergence  to  steady  state.  A  steady-state  solution  was  possible  because 
the  unsteady  shedding  phenomena  associated  with  bluff  body  flow  were  substantially 
damped  by  the  Spalart-Allmaras  RANS  turbulence  model. 


(3.3) 
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3.5.2  Detached-Eddy  Simulation  Specific  Setup.  Once  a  steady-state  so¬ 
lution  was  attained  using  strictly  RANS  methods,  the  DES  turbulence  model  was 
selected  on  the  Viscous  panel.  This  action  automatically  set  the  solver  to  an  un¬ 
steady,  second-order  temporal  discretization  employing  a  bounded,  central- difference 
discretization  of  the  momentum  equation.  The  default  DES  constant  of  Cdes  —  0.65 
suggested  by  Spalart  was  accepted.  The  Spalart-Allmaras  turbulence  model  is  used 
exclusively  as  the  RANS  model  for  a  DES  in  Fluent™  6.2.16. 

Finally,  when  using  Fluent’s™  implicit  solver,  which  employs  dual  time-stepping, 
the  maximum  number  of  sub-iterations  for  each  physical  time  step  must  be  deter¬ 
mined.  Essentially,  each  time  step  is  a  steady-state  computation  and  should  be  suffi¬ 
ciently  converged  before  incrementing  the  physical  time  step.  A  sub-iteration  conver¬ 
gence  study  using  the  9.9  million  cell  grid  was  conducted  to  determine  an  appropriate 
maximum  number  of  sub-iterations  for  the  solver  during  DES.  Figure  3.4  shows  resid¬ 
ual  plots  over  several  timesteps.  The  residuals  begin  to  show  asymptotic  behavior 
after  50  sub-iterations.  However,  for  100  sub-iterations  per  timestep,  the  residuals 
were  completely  asymptotic.  Thus,  during  data  collection  runs,  70  sub-iterations  was 
deemed  more  than  sufficient. 

3.6  Visualization  Tool  Kit  and  Post-Processing 

Large  grid  sizes  and  small  timesteps  lead  to  significant  data  storage  and  post¬ 
processing  burdens.  Fluent™  Data  hies  for  a  single  time  step  contain  the  how  proper¬ 
ties  at  at  every  point  in  the  grid  and  ballooned  in  size  to  2.9  Gigabytes  for  the  14  mil¬ 
lion  cell  grid.  Using  Fluent’s™  built-in  data  output  format,  collecting  four  seconds  of 
time-accurate  data  for  a  single  case  would  require  approximately  1.1  Terabytes.  An 
alternative  to  Fluent  is™  data  storage  method  was  the  Fluent  VTK  Extractor  (FVE) 
which  directly  accesses  data  from  the  Fluent™  arrays  at  runtime  and  outputs  user- 
specified  geometry  and  data  to  a  Visualization  Tool  Kit  (VTK)  XML  hie.  [15]  VTK 
is  an  open  source,  3-D  graphics,  image,  and  visualization  software  application  which 
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Figure  3.4:  Residual  Output,  9.9M  cell  grid,  1  m/s  inlet  velocity  (a)  50  Sub¬ 

iterations  per  timestep  (b)  100  Sub-iterations  per  timestep 
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Figure  3.5:  VTK  bounding  box;  only  data  inside  the  box  is  stored  at  runtime; 

boundary  layer  and  Strouhal  number  were  computed  for  the  numbered  cylinders 

supports  scalar,  vector,  and  tensor  operations,  as  well  as  mesh  smoothing  and  cutting 
functions  for  filtering  data.  [20] 

The  FVE  is  implemented  as  a  Fluent™  User-Defined,  Compiled  Function.  The 
user  first  creates  C  functions  containing  instructions  for  specifying  bounding  boxes, 
cut-planes,  and  point  probes  used  to  sample  scalars  of  interest  such  as  pressure, 
density,  and  velocity  at  runtime.  The  C  hie  and  compiled  FVE  libraries  are  selected 
on  the  Fluent™  Compiled  Functions  menu  and  then  subsequently  built  and  loaded. 
Finally,  at  runtime,  Fluent™  calls  the  FVE  user-defined  function  and  outputs  a  VTK 
format  data  hie  for  each  time  step.  This  method  saved  tremendous  storage  space 
since  only  regions  of  interest  in  the  grid  are  stored,  and  the  number  of  scalars  in  the 
output  is  reduced  to  those  selected  by  the  user.  Data  storage  savings  of  up  to  99% 
and  runtime  reductions  of  up  to  15%  are  possible  when  using  FVE.  [15]  Figure  3.5 
illustrates  the  size  and  location  of  the  data  subset  saved  at  runtime  using  FVE. 

Because  the  fluid  is  an  incompressible,  ideal  gas  formulation  employing  the 
Segregated  Solver,  density  is  treated  as  a  constant.  To  minimize  output,  only  static 
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pressure  and  x,  y,  and  z  velocity  were  chosen  for  output  during  DES.  With  these  four 
scalars,  virtually  any  other  flow  quantity  can  be  computed  during  post-processing. 

Because  the  simulation  was  run  in  parallel,  the  FVE  utility  wrote  an  output  hie 
for  each  partition  of  the  grid  at  each  time  step.  In  order  to  view  the  entire  dataset,  the 
information  must  be  stitched  back  together  using  a  locally  developed  utility  written  in 
Python  known  as  the  Fuser.  The  Fuser  uses  VTK  function  calls  to  identify  the  shared 
nodes  of  each  partition  and  output  a  single,  unified  grid  and  data  set  for  visualization 
and  data  reduction. 

Visualization  of  the  VTK  output  was  accomplished  with  MayaVi,  an  open- 
source  visualization  program  written  in  Python.  Using  MayaVi,  each  of  the  stored 
scalars  can  be  viewed  and  interpreted  with  a  variety  of  techniques  including  cut-planes 
and  point  probes.  [16]  MayaVi  is  limited,  however,  as  far  as  computing  advanced  how 
dynamics  such  as  vorticity  and  path  lines  in  place.  For  this  reason,  a  Python  script, 
vtk2ensight,  was  developed  to  convert  the  VTK  hies  to  Ensight™  format.  Vtk2ensight 
uses  VTK  5.0  to  write  each  of  the  scalar  arrays  in  the  VTK  output  to  time-accurate 
Ensight™  data  hies. 

Ensight™,  a  commercial  engineering  and  scientific  visualization  application,  is 
able  to  recognize  transient  data  sets  and  import  each  time  step  into  a  single  session. 
This  allows  the  user  to  examine  changes  in  the  data  over  time  using  animation  and 
particle  path-lines.  Furthermore,  Ensight™  has  the  ability  to  readily  compute  and 
visualize  vorticity,  velocity  vectors,  and  other  derived  how  characteristics.  [6] 

Animations  are  an  extremely  powerful  tool  for  garnering  an  understanding  of 
unsteady  how  dynamics.  Thus,  animated  visualizations  of  derived  variables  such  as  ve¬ 
locity  magnitude,  velocity  vectors,  and  vorticity  were  constructed  using  Ensight ’s™  vari¬ 
able  calculator  and  transient  data  hip-books.  Path-lines,  the  locus  of  points  describing 
the  movement  of  a  massless  particles  through  an  unsteady  velocity  held  are  another 
tool  available  inside  Ensight.  Naturally,  a  visualization  capability  such  as  this  was 
extraordinarily  helpful  in  ascertaining  mechanisms  for  droplet  transport. 
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Examination  of  instantaneous  data  in  the  form  of  animations  and  path-lines 
is  a  valuable  tool,  however,  time-averaged  data  considers  long-period  bulk  motion. 
Computing  time-averaged  data  required  the  development  of  another  Python  script 
with  VTK  function  calls.  The  time-averaging  utility  read  in  each  scalar  array  from 
VTK  XML  data  hies,  summed  the  data,  and  divided  by  the  number  of  timesteps 
before  finally  outputting  the  data  to  a  new  VTK  XML  hie. 

Determining  spectral  content  in  a  how  at  a  point  is  the  hrst  step  in  computing 
the  Strouhal  number  (Section  2.1).  Extracting  data  collected  via  a  point  probe  over 
time  necessitated  the  development  of  another  Python  script.  Scalar  data  from  time- 
accurate  VTK  data  hies  was  interpolated  from  cell-centered  quantities  to  a  point 
requested  by  the  user  and  written  to  an  ASCII  hie.  The  resultant  record  of  time- 
accurate  data  at  a  point  in  the  how,  formatted  as  space-delimited  ASCII,  was  then 
readily  accessible  to  Matlab™  for  performing  spectral  analysis  and  computing  tur¬ 
bulence  statistics.  Power  Spectral  Density  (PSD)  plots  of  the  data  were  created  to 
compute  the  Strouhal  number  using  Equation  (2.1).  Turbulence  statistics  such  as 
mean  velocity,  RMS  velocity,  and  turbulence  intensity  were  also  checked. 

To  conclude,  the  synthesis  of  time-accurate  methods  such  as  animated  visualiza¬ 
tions,  path-lines,  and  spectral  content  combined  with  time-averaged  data  yielded  the 
necessary  tools  to  reduce  and  characterize  the  turbulent  how-held  inside  the  cylinder 
array  and  search  for  explanations  of  droplet  transport. 
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IV.  Results 


The  results  presented  in  this  chapter  summarize  the  outcome  of  RANS  and 
Detached-Eddy  Simulations  from  the  standpoint  of  model  validation  and  droplet 
transport  mechanisms  related  to  fire  suppression.  In  the  vein  of  model  validation,  the 
numerical  simulations  focused  on  boundary  layer  and  vortex  shedding  physics  for  a 
single  cylinder  and  an  array  of  cylinders.  Examination  of  the  data  from  the  perspec¬ 
tive  of  droplet  transport  was  conducted  exclusively  on  the  fine  cylinder  array  grid. 

4-1  Single  Cylinder  DES  Validation 

4-1.1  Boundary  Layer  Resolution.  Data  gathered  from  both  RANS  and 
DES  methods  supported  the  verification  of  the  boundary  layer  profile.  Wall  y+,  a 
non-dimensional  parameter  used  to  quantify  the  quality  of  near-wall  grid  spacing, 
was  plotted  for  each  case  to  confirm  Spalart-Allmaras  performance  in  the  boundary 
layer.  When  using  the  Spalart-Allmaras  equation  to  model  the  boundary  layer  and 
avoid  wall-functions,  the  Fluent™  manual  recommends  y+  values  on  the  order  of  one 
or  less  at  wall  boundaries.  [9]  Figure  4.1  gives  RANS  and  DES  wall  y+  values  for  the 
single  cylinder  grid  at  =  1  m/s.  Notably,  y+  <  1  was  satisfied  everywhere  along 
the  surface  of  the  cylinder.  Similarly,  y+  trends  for  the  5  m/s  free-stream  case  were 
also  observed. 

In  addition  to  inspecting  y+,  the  curvature  of  the  boundary  layer  profile  was  also 
qualitatively  analyzed  to  assess  the  quality  of  the  boundary  layer  model.  Turbulent 
boundary  layers  along  a  line  described  by  Figure  4.3  are  clearly  depicted  in  Figure  4.2, 
demonstrating  adequate  grid  spacing.  The  results  of  the  wall  y+  variables  coupled 
with  the  boundary  layer  profiles  offers  sufficient  evidence  to  accept  the  boundary  layer 
as  well-resolved  and  sufficient  for  computing  an  accurate  numerical  solution. 

4-1.2  Spectral  Content.  Frequency  data  collected  at  a  point  described  by 
Figure  4.4  was  used  to  compute  the  Strouhal  number  for  the  single  cylinder  2.2  mil¬ 
lion  cell  grid.  The  power  spectral  density  (PSD)  plot  for  the  single  cylinder  case, 
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Wall  y+  for  a  Single  Cylinder 
1  m/s 


Angle  (deg) 
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Figure  4.1:  RANS  and  DES  comparison  of  wall  y+  values  for  a  single  cylinder, 

U0 0  =  1  m/s,  2.2  M  cells 


Boundary  Layer  Profile  at  90°  for  a  Single  Cylinder 
1  m/s 


Figure  4.2:  RANS  and  DES  boundary  layer  profile  comparison  for  a  single  cylinder 
at  90°,  Uoo  =  1  m/s,  2.2  M  cells;  Figure  4.3  describes  the  measurement  location 
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Figure  4.3:  Boundary  layer  profile  measurement  location;  the  arrow  describes  the 
direction  of  increasing  ynorm 


Figure  4.4:  Frequency  measurement  location  for  computing  Strouhal  number 
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Power  Spectrum  for  a  Single  Cylinder  at  Point  1  (X  Velocity) 
1  m/s 


(a) 


Power  Spectrum  for  Single  Cylinder  at  Point  1  (X  Velocity) 
5  m/s 


(b) 

Figure  4.5:  Power  spectrum,  single  cylinder,  2.2  M  cells  (a)  =  1  m/s  (b)  Uoo  = 

5  m/s;  Figure  4.4  describes  the  measurement  location 
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Wall  y+  for  Cylinder  Array  (9.9M  Cells) 
5  m/s 


Angle  (deg) 

Figure  4.6:  RANS  and  DES  comparison  of  wall  y+  values,  cylinder  array,  UQ 0  = 
5  m/s,  9.9  M  cells 

Figure  4.5  (a),  shows  a  strong  peak  at  8  Hz  corresponding  to  a  Strouhal  number  of 
St  =  0.203  (Re  =  1,575)  matching  nicely  with  the  accepted  value  of  St  =  0.2  for 
100  <  Re  <  105.  The  higher  speed  case,  Uoo  =  5  m/s,  depicted  in  Figure  4.5  (b)  pro¬ 
vides  no  obvious  shedding  frequency  indicating  the  grid  was  too  coarse  to  adequately 
capture  the  higher  frequency  mode  associated  with  a  faster  free-stream  velocity.  While 
the  single  cylinder  grid  was  not  refined  to  resolve  the  shedding  in  the  high  speed  case, 
the  results  influenced  the  refinement  of  the  cylinder  array  grid  to  ensure  shedding  was 
captured. 

4-2  Exploration  of  the  Circular  Cylinder  Array 

4-2.1  Boundary  Layer  Resolution.  Investigation  of  boundary  layer  behavior 
on  the  cylinder  array  grids  was  accomplished  in  the  same  manner  as  the  analysis 
performed  on  the  single  cylinder  grid.  A  plot  of  y+  values  for  the  cylinder  array  is 
given  in  Figure  4.6.  Although  the  figure  is  crowded,  the  the  maximum  wall  y+  values 
are  the  key  feature  exhibited  by  the  plot.  Once  more,  the  Spalart-Allmaras  criteria 
of  wall  y+  <  1  was  met  throughout  the  array. 
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Figure  4.7:  VTK  bounding  box  illustrating  the  designation  of  certain  cylinders  in 
the  array 

Figure  4.8  qualitatively  compares  the  boundary  layer  profiles  of  Cylinders  1, 
2,  and  3  designated  in  Figure  4.7,  taken  along  a  line  described  by  Figure  4.3.  Like 
the  single  cylinder  case,  the  boundary  layer  velocity  profiles  at  U0 0  =  5  m/s  on  the 
9.9  million  cell  array  grid  presented  well-captured  turbulent  velocity  profiles.  Since 
the  variation  of  free-stream  velocity  and  grid  density  did  not  appreciably  influence 
the  behavior  of  the  boundary  layer,  only  a  single  case  was  illustrated  in  Figure  4.8. 

4-2.2  Spectral  Content.  Even  though  the  Strouhal  number  corresponding  to 
a  circular  cylinder  in  an  array  is  not  correlated  to  Reynolds  number  in  the  same  way 
as  it  is  for  a  single  cylinder,  the  frequency  content  and  Strouhal  number  for  Cylinder  1 
were  nonetheless  computed  at  a  point  indicated  in  Figure  4.4.  Cylinder  1  was  chosen 
because  the  oncoming  flow  is  free-stream,  unlike  the  flow  seen  by  Cylinders  2  and  3, 
which  lie  in  the  unsteady,  turbulent  wakes  of  the  upstream  cylinders  and  experience 
mean  local  flow  velocities  as  high  as  twice  the  free-stream.  The  PSD  plots  of  Figure  4.9 
from  the  9.9  million  cell  grid  revealed  a  faint  peak  near  6  Hz  for  =  1  m/s  and 
no  discernible  shedding  frequency  for  =  5  m/s.  The  lack  of  a  shedding  frequency 
for  the  5  m/s  free-stream  case  indicates  that  the  9.9  million  cell  grid  was  too  coarse 
to  detect  the  flow  physics  of  interest.  Nevertheless,  the  enhanced  spatial  resolution 
of  the  14  million  cell  grid  satisfactorily  captured  the  Karman  vortex  shedding  which 
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Boundary  Layer  Profile  at  90°  for  Cylinder  1  (9.9M  Cells) 
5  m/s 


(a) 


Boundary  Layer  Profile  at  90°  for  Cylinder  2  (9.9M  Cells)  Boundary  Layer  Profile  at  90°  for  Cylinder  3  (9.9M  Cells) 

5  m/s  5  m/s 


(b)  (c) 

Figure  4.8:  RANS  and  DES  cylinder  array  boundary  layer  profile  comparison  at  90°, 
Uoo  =  5  m/s,  9.9  M  cells;  Figure  4.3  describes  the  measurement  location;  (a)  Cylin¬ 
der  1  (b)  Cylinder  2  (c)  Cylinder  3 
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Figure  4.9: 

1  m/s  (b)  U, 


Power  Spectrum  for  Cylinder  Array  at  Point  1  (X  Velocity) 
1  m/s 


(a) 


Power  Spectrum  for  Cylinder  Array  at  Point  1  (X  Velocity) 
5  m/s 

101  - . - . - 1 - 1 - 1 - , - . - 


Frequency  (Hz) 


(b) 

Power  spectrum,  cylinder  array  at  cylinder  1,  9.9  M  cells  (a)  U 
=  5  m/s;  Figure  4.4  describes  the  measurement  location 
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Figure  4.10: 
1  m/s  (b)  Uc 


Power  Spectrum  for  Cylinder  Array  at  Point  1  (X  Velocity) 
1  m/s 


(a) 


Power  Spectrum  for  Cylinder  Array  at  Point  1  (X  Velocity) 
5  m/s 


(b) 

Power  spectrum,  cylinder  array  at  cylinder  1,  14  M  cells  (a) 
=  5  m/s;  Figure  4.4  describes  the  measurement  location 
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is  evident  in  the  power  spectrum  plots  of  Figure  4.10.  At  Uqo  =  1  m/s,  a  10  Hz 
signal  dominates  the  flow,  and  at  Uoo  =  5  m/s,  a  20  Hz  peak  represents  the  shedding 
frequency.  Strouhal  numbers  based  on  the  14  million  cell  grid  data  were  computed 
at  Uoo  =  1  m/s  and  U0 0  =  5  m/s  to  be  St  =  0.25  and  St  =  0.11  respectively.  The 
deviation  from  the  accepted  St  =  0.2  for  a  single  cylinder  is  likely  due  to  influence  of 
adjacent  cylinders  shedding  asynchronously  with  Cylinder  1. 

The  results  of  the  power  spectral  density  plots  of  the  single  cylinder  and  the 
cylinder  array  demonstrated  that  the  objective  cell  spacing  of  A0  =  0.2  in  the  LES 
region  was  insufficient  to  describe  the  primary  flow  physics  of  interest.  Although 
the  14  million  cell  mesh  was  intended  to  investigate  solution  grid-independence,  the 
simulation  remained  a  success  because  the  14  million  cell  grid  enabled  the  study  of 
droplet  transport  mechanisms. 

4-3  Droplet  Transport  Mechanisms 

Initially,  the  search  for  droplet  transport  and  trapping  mechanisms  focused  on 
span-wise  flow  and  vortices  capable  of  entraining  particles.  The  first  hypothesis  put 
forth  in  Section  1.3  contended  that  non-uniform,  span-wise  shedding  generated  waves 
pushing  entrained  droplets  toward  the  side  walls  of  the  tunnel.  After  deposition  on 
the  walls,  the  droplets  would  drip  and  pool  on  the  tunnel  floor. 

Strong,  non-uniform  shedding  along  the  span  of  a  single  cylinder  is  typically 
manifested  by  the  pinching  and  folding  of  vortex  cores  in  the  wake  illustrated  by 
Figure  4.11.  For  the  single  cylinder  case,  weak,  non-uniform  shedding  was  observed 
and  is  depicted  in  Figure  4.12  and  Figure  4.13.  As  time  progresses,  the  vortex  core 
represented  by  (a)  and  (b)  in  Figure  4.12  bends  and  separates  into  two  distinct  vortex 
cores  at  (c)  and  (d)  causing  weak  span-wise  flow  to  develop  briefly.  The  y-velocity 
contours  shown  in  Figure  4.13  correspond  to  the  the  vortex  core  evolution  in  Fig¬ 
ure  4.12.  The  y-velocity  plots,  taken  0.25H  downstream  of  the  cylinder,  depict  weak 
asynchronous  span- wise  shedding.  As  (a)  and  (b)  of  Figure  4.13  are  trending  down, 
(d)  and  (e)  indicate  the  opposite.  While  the  non-uniform  shedding  and  span-wise 
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Vortex  core 


Figure  4.11:  Expected  vortex  core  time  evolution  during  strong,  non-uniform,  span- 
wise  shedding 


vortic-ity  of  the  single  cylinder  case  was  detectable,  the  span-wise  flow  was  quite  weak 
and  not  capable  of  supporting  span-wise  droplet  transport. 

The  same  analysis  of  vortex  cores  and  y-velocity  distribution  on  the  cylinder 
array  grid  yielded  no  evidence  whatsoever  of  span-wise  vorticity  or  non-uniform  shed¬ 
ding.  Figure  4.14  shows  a  vortex  core  computation  representative  of  the  cylinder 
array  which  is  remarkably  more  uniform  and  stratified  than  that  single  cylinder  case. 
Y-velocity  plots  on  a  cut-plane  0.25 D  downstream  of  the  second  row  of  cylinders 
(Figure  4.15)  further  reinforces  the  absence  of  non-uniform  shedding  and  the  lack  of 
span-wise  flow  inside  the  array.  In  fact,  there  is  almost  no  change  in  the  y-velocity  dis¬ 
tribution  along  the  span  of  the  cylinders.  The  close  proximity  of  the  cylinders  seems 
to  force  the  vortex  cores  to  remain  aligned  with  the  cylinder  span.  Thus,  span-wise 
droplet  transport  toward  the  tunnel  walls  was  not  witnessed  using  this  model  for  the 
conditions  specified. 


4-3.1  Particle  Transit  Time  Analysis.  With  the  span-wise  motion  hypoth¬ 
esis  disproved,  the  second  hypothesis  described  in  Section  1.3  was  explored.  A  first- 
order  inertial  estimate,  originally  discussed  in  Section  2.4.1,  was  conducted  using  data 
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Figure  4.12:  Vortex  core  time  evolution  during  weak,  non-uniform,  span-wise  shed¬ 
ding,  single  cylinder,  Uoo  =  1  m/s 
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Figure  4.13:  Y-Velocity  visualization,  Uoo  =  1  m/s,  single  cylinder,  x  =  0.25D 
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Figure  4.14:  Instantaneous  vortex  core  visualization,  (7^  =  1  m/s,  cylinder  array, 
14  M  cells 


from  the  14  million  cell  grid  Detached-Eddy  Simulations.  In  the  most  conservative 
case,  a  droplet  could  fall  as  far  as  1  inch  (0.0254  m ),  the  vertical  distance  between 
cylinders  in  the  same  row,  before  impacting  the  clutter.  For  the  1  m/s  free-stream 
case,  the  average  carrier  fluid  velocity  in  the  spaces  between  the  cylinders  of  the  array 
was  approximately  1.75  m/ s.  Applying  only  the  acceleration  of  gravity  to  the  droplet, 
the  time  required  for  a  droplet  to  fall  and  impact  another  cylinder  was  t  =  0.072  s. 
However,  at  a  mean  velocity  of  1.75  m/s  the  time  required  for  a  particle  to  traverse 
7  inches  (0.1778  m),  the  horizontal  dimension  of  the  array,  was  t  =  0.10  s.  Thus, 
the  likelihood  of  a  droplet  impacting  one  of  the  cylinders  before  exiting  the  array 
was  quite  high  in  the  low-speed  case.  In  contrast,  the  average  velocity  between  the 
cylinders  for  the  5  m/s  free-stream  case  was  8.5  m/s.  A  droplet  entrained  in  this  flow 
could  negotiate  the  array  in  approximately  t  =  0.02  s,  which  is  3.6  times  faster  than 
the  time  required  for  the  same  particle  to  fall  into  another  cylinder.  Although  this 
analysis  is  quite  simple,  it  indicates  that  velocity  and  inertia  play  the  dominating  role 
in  governing  whether  or  not  an  entrained  particle  impacts  a  cylinder  or  is  whisked 
through  the  clutter. 
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Figure  4.15:  Y-Velocity  visualization,  Uoo  =  1  m/s,  cylinder  array,  x  =  0.25D, 

14  M  cells 
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Figure  4.16:  Streamline  and  Stokes  number  visualization  at  simulation  t  =  3.77  s, 
Uoo  =  1  m/s,  cylinder  array,  14  M  cells 


4-3.2  Stokes  Number  Analysis.  Investigating  whether  or  not  the  droplets 
were  likely  to  impact  the  clutter  based  on  unsteady  phenomena  as  well  as  the  time- 
averaged  flow  characteristics  was  accomplished  by  composing  animations  of  stream¬ 
lines  superimposed  on  contour  plots  of  Stokes  number.  This  was  inspired  by  the 
discussion  of  Section  2.4.2  which  described  the  effect  of  Stokes  number  on  entrained 
particle  motion.  The  Stokes  number  computation  assumed  water  droplet  particles  of 
approximately  90  fim  in  diameter,  a  likely  size  for  droplets  emanating  from  a  spray 
nozzle.  [9] 

Figure  4.16  contains  a  representative  instantaneous  Stokes  number  contour  plot 
with  streamlines  superimposed  for  the  1  m/s  free-stream  case.  Overall,  the  Stokes 
number  throughout  the  flow-field  depicted  in  Figure  4.16  is  relatively  low.  Much  of 
the  flow  between  the  cylinders  is  colored  green  which  corresponds  to  Stokes  numbers 
on  the  order  of  one;  also  the  peak  Stokes  number  of  the  flow,  Stk  =  2.6,  is  relatively 
low.  The  vortex  shedding  of  the  cylinders  in  Figure  4.16  induced  significant  streamline 
curvature  in  the  wake  where  the  corresponding  Stokes  number  is  on  the  order  of  one  or 
less  indicating  that  any  dispersed  particles  would  likely  track  along  the  streamlines. 
Furthermore,  streamlines  in  the  low-speed  case  illustrated  by  Figure  4.16  (a)  pass 
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Figure  4.17: 
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Streamline  and  Stokes  number  evolution,  Uoo  =  1  m/s ,  cylinder  array, 
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Figure  4.18:  Streamline  and  Stokes  number  visualization  at  simulation  t  =  3.75  s, 
U0 o  =  5  m/s,  cylinder  array,  14  M  cells 

through  or  very  near  shed  vortices  where  the  Stokes  number  is  small.  Thus,  droplets 
following  these  streamlines  could  easily  be  entrained  inside  the  vortex  and  carried 
along.  Animations  of  the  Stokes  number,  summarized  by  Figure  4.17,  clearly  show 
that  vortices  shed  from  the  second  row  of  cylinders  in  the  1  m/s  case  often  directly 
impact  the  final  row  of  cylinders.  Figure  4.17  (a)  points  out  a  vortex  shed  from 
Cylinder  2  which  impacted  Cylinder  3  (Figure  4.7  for  cylinder  numbers).  Naturally, 
any  droplets  entrained  in  a  vortex  that  impacts  a  cylinder  would  likely  be  deposited 
on  the  cylinder,  accumulate,  and  drip  to  the  floor. 

Conversely,  in  the  5  m/s  free-stream  case,  shown  in  Figures  4.18  and  4.19, 
there  are  fewer  highly  curved  streamlines  in  the  wake,  and  those  that  are  present  lie 
in  regions  of  Stokes  numbers  greater  than  one  (Figure  4.19  (a)),  meaning  droplets 
attempting  to  follow  the  carrier  fluid  would  fail  to  negotiate  the  streamline  curvature 
and  instead  follow  a  more  ballistic  path.  Also,  in  the  5  m/s  free-stream  case,  the 
animations  showed  far  fewer  vortices  shed  from  the  second  row  of  cylinders  impacting 
the  third  row.  Thus  any  droplets  that  were  under  the  influence  of  a  Karman  vortex 
would  not  likely  be  deposited  on  a  cylinder.  As  a  result,  fewer  droplets  in  the  high 
speed  case  would  be  influenced  by  the  unsteady  shedding  phenomenon  and  therefore 
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Increasing  Time 


Figure  4.19:  Streamline  and  Stokes  number  evolution,  =  5  m/s,  cylinder  array, 
14  M  cells 
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the  motion  of  entrained  droplets  will  likely  be  governed  primarily  by  average  flow- 
field  characteristics.  This,  coupled  with  the  transit  time  analysis  of  Section  4.3.1, 
supports  the  observations  of  by  Disimile  et  al.  [7]  and  serves  as  an  explanation  for  the 
entrapment  and  transport  of  droplets  inside  idealized  engine  nacelle  clutter. 
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V.  Conclusions 


Many  of  the  subtleties  of  CFD  and  turbulence  modeling  and  analysis  were 
discovered  during  the  process  of  investigating  droplet  transport  mechanisms 
through  an  array  of  circular  cylinders. 

5.1  Detached-Eddy  Simulation 

The  Detached-Eddy  Simulation  turbulence  modeling  technique  was  successfully 
demonstrated  as  an  effective  tool  for  computing  large  turbulent  flow  structures  in  low- 
Reynolds  number  flows  at  a  relatively  modest  computational  expense  when  compared 
with  LES.  In  cases  such  as  the  cylinder  array,  where  the  unsteady  phenomena  of 
interest  are  large  eddies  outside  of  the  boundary  layer,  DES  is  an  excellent  numerical 
method. 

The  most  obvious  CFD  lesson  in  this  research  was  stated  best  by  Spalart,  “...any 
unsatisfactory  result  reported  to  the  author  [Spalart]  is  due  to  the  user’s  failure  to  run 
on  a  fine  enough  grid.”  [21]  Without  a  doubt,  failure  to  use  a  sufficiently  refined  grid 
and  appropriately  small  timesteps  resulted  in  poorly  resolved  transient  behavior  on 
the  9.9  million  cell  grid.  Fortunately,  the  14  million  cell  grid,  originally  only  intended 
as  a  grid  convergence  test,  adequately  captured  the  bulk  shedding  phenomenon  of  the 
array  for  both  inlet  velocity  conditions.  In  the  future,  more  care  should  be  given  to 
spatial/temporal  error  balancing.  The  guidance  Spalart  gives  for  estimating  the  time 
step,  At  =  A0 /Umax,  discussed  in  section  2.3.4,  should  be  treated  only  as  an  initial 
estimate.  The  user  should  experiment  freely  with  space-time  balancing  and  err  on 
the  side  of  small  time  steps. 

5.2  Droplet  Transport 

Two  hypotheses  were  postulated  at  the  outset  of  this  research  to  explain  the 
droplet  transport  and  entrapment  phenomena  observed  by  Disimile  et  al.  [7].  The  first 
hypothesis,  which  supposed  span-wise  flow  developed  in  the  cylinder  wake  region  due 
to  non-uniform  cylinder  shedding,  was  not  observed  using  this  model.  Instead,  this 
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research  successfully  demonstrated  the  veracity  of  a  second  hypothesis,  which  posited 
the  dominance  of  the  velocity-inertia  relationship  as  the  governing  mechanism  for 
dispersed-phase  particles  negotiating  obstacles. 

For  slower  average  velocities  inside  the  array,  entrained  suppressant  droplets 
will  likely  fall  or  be  carried  into  neighboring  cylinders  before  successfully  passing 
through  the  cylinder  array.  Conversely,  particles  entrained  in  faster  flow-fields  are 
able  to  transit  the  array  before  the  time  required  to  fall  into  an  obstacle  has  elapsed. 
Furthermore,  low  Stokes  number  regions  capable  of  capturing  entrained  droplets  were 
revealed  by  animations  of  vortex  shedding  in  the  cylinder  wake.  Vortices  that  were 
capable  of  entraining  droplets  after  shedding  from  the  second  row  of  cylinders  often 
directly  collided  with  the  final  row  of  cylinders.  Such  a  collision  would  certainly 
deposit  the  droplets  carried  by  the  vortex  on  the  impacted  cylinder  where  the  droplets 
would  accumulate  and  eventually  fall  to  the  tunnel  floor  as  super-droplets  too  large 
for  entrainment. 

From  the  perspective  of  transporting  fire  suppressant  past  bluff  body  clutter 
to  a  downstream  fire,  high-airspeed  flows  will  clearly  be  more  effective  for  entraining 
and  transporting  suppressant  droplets  through  clutter  without  impact.  However  if 
the  cylinders  of  this  research  represent  fuel  lines,  the  goal  of  suppressant  transport 
would  not  be  swift  conveyance  through  the  array.  If  one  of  the  fuel  lines  in  the  array 
burst  and  began  spraying  fuel  into  the  cylinder  wake,  a  bluff-body  stabilized  turbu¬ 
lent  diffusion  flame  will  form  if  the  fuel  were  ignited.  [1]  Under  these  circumstances, 
particles  dispersed  from  a  suppressant  spray  nozzle  into  a  high-speed  co-flow  will  not 
likely  reach  the  cylinder’s  recirculation  region,  resulting  in  a  fire  that  is  very  difficult 
if  not  impossible  to  extinguish. 

5.3  Future  Research 

Future  research  regarding  the  cylinder  array  should  attempt  to  establish  the 
effect  of  free-stream  turbulence  intensity  on  the  flow  held  inside  the  cylinder  array 
by  using  a  DES  based  on  a  R.ANS  model  other  than  Spalart-Allmaras  so  that  free- 
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stream  disturbances  can  propagate  towards  the  clutter  instead  of  being  damped  to 
zero.  Another  variable  of  interest,  which  was  investigated  by  Disimile  et  al.  [7],  was 
the  horizontal  cylinder  spacing  between  the  rows  of  cylinders.  Disimile  et  al.  logically 
found  that  tighter  cylinder  spacing  equated  to  less  droplet  transport.  However,  further 
investigation  of  the  spacing  variable  might  provide  additional  insight  into  the  overall 
droplet  transport  problem.  Finally,  a  simulation  employing  Fluent’s™  dispersed- 
phase  model  should  be  attempted  to  individually  track  the  particles  of  the  entrained 
suppressant  droplets. 

Fluent™  is  a  nominally  second-order  solver,  meaning  that  the  order  of  the  trun¬ 
cation  error  of  the  discretization  of  the  governing  equations  is  second-order.  Therefore, 
as  the  grid  is  refined  and  a  new  solution  is  computed,  the  truncation  error  should  de¬ 
crease  by  a  power  of  two.  Logically,  it  follows  that  increased  grid  refinement  is  required 
to  reach  an  acceptably  accurate  solution  for  a  second-order  solver  than  a  higher-order 
solver.  As  a  result,  the  requirement  of  finer  meshes  for  DES  in  Fluent™  increases 
both  the  computational  and  data  storage  requirements  to  support  the  solution.  Of 
course,  using  a  solver  with  an  order  of  accuracy  higher  than  second-order  will  likely 
require  using  a  structured  grid  to  facilitate  the  larger  discretization  stencil.  Introduc¬ 
ing  structured  grid  topology  to  the  problem  ushers  in  a  whole  new  set  of  obstacles  to 
DES  grid  generation,  however. 

The  initial  condition  for  a  Detached-Eddy  Simulation  is  usually  a  steady-state 
RANS  solution.  Naturally  once  the  DES  begins,  transient  artifacts  resulting  from  the 
RANS  to  DES  transition  will  persist  in  the  solution  and  distort  the  data.  However, 
the  length  of  time  associated  with  the  transients  varied  more  than  anticipated  from 
one  solution  to  the  next.  As  a  result,  computational  effort  (time)  was  unnecessarily 
expended  computing  and  storing  data  which  was  ultimately  discarded.  In  the  future, 
several  seconds  of  unsteady  data  should  be  computed  without  utilizing  the  FVE  to 
accelerate  obtaining  a  fully-developed  unsteady  flow-held. 
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Finally,  demonstrating  grid  independence  of  the  solution  was  not  accomplished 
in  this  study  because  the  fine  mesh  intended  for  this  purpose  was  actually  the  min¬ 
imum  grid  resolution  necessary  to  capture  the  pertinent  flow  physics.  A  25  million 
cell  grid  was  constructed  to  meet  the  needs  of  the  grid  convergence  study,  but  the 
grid  could  not  be  exported  to  Fluent™  because  Gridgen ™,  a  32-bit  application,  was 
unable  to  allocate  enough  memory.  If  the  25  million  cell  grid  can  be  exported,  the 
future  researchers  should  be  cognoscente  of  the  associate  runtime  for  such  a  large 
grid.  The  average  runtime  for  the  14  million  cell  grid  was  188  hours  or  7.9  days  to 
collect  two  seconds  of  data  when  employing  the  Fluent  VTK  Extractor  (FVE).  Based 
on  this,  the  runtime  for  the  25  million  cell  grid  will  be  approximately  two  weeks  to 
compute  two  seconds  of  data. 

5.4  Summary 

In  conclusion,  this  research  modeled  the  idealized  engine  nacelle  clutter  of  Dis¬ 
imile  et  al.  in  an  effort  to  ascertain  flow-held  dynamics  capable  of  explaining  the 
failure  of  water  droplets  to  transit  an  array  of  cylinders  in  low-speed  co-how.  Numer¬ 
ical  simulations  employing  Detached-Eddy  Simulation  were  conducted  to  resolve  the 
turbulent  cylinder  wakes  and  examine  unsteady  how-held  dynamics  inside  the  array. 
Ultimately,  this  study  determined  that  droplets  emanating  from  a  spray  nozzle  in 
co-how  would  follow  paths  governed  by  the  relationship  between  mean  carrier  fluid 
velocities,  droplet  inertia,  and  local,  unsteady,  low-Stokes  number  regions.  For  low 
free-stream  velocities  (~  1  m/s)  suppressant  droplets  will  likely  become  entrained  in 
shed  vortices  as  the  droplets  pass  through  low-Stokes  number  portions  of  the  how.  In 
the  low-speed  case,  vortices  shed  from  the  second  row  of  cylinders  typically  impacted 
the  third  row  of  cylinders  directly;  this  action  would  likely  deposit  any  entrained 
droplets  onto  the  cylinder  where  they  would  aggregate  and  drip.  On  the  other  hand, 
droplets  entrained  in  the  higher-speed  how  (5  m/s)  would  not  be  effected  by  the 
carrier  fluid's  streamlines  because  the  Stokes  number  in  between  the  cylinders  was 


60 


substantially  higher.  Thus,  these  particles  would  likely  transit  the  array  with  ease  at 
mean  carrier  fluid  velocities. 
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Appendix  A.  Ancillary  CFD  Discussions 

A  few  subtleties  of  CFD  are  discussed  in  this  section  outside  of  the  main  document 
because  they  do  not  immediately  relate  to  the  overall  objective:  determining  mecha¬ 
nisms  of  droplet  transport  inside  engine  nacelles.  However,  the  information  presented 
here  was  discovered  during  the  course  of  the  research  effort. 

A.l  Modified  Turbulent  Viscosity  and  Turbulence  Intensity 

A  DES  based  on  the  Spalart-Allmaras  model  cannot  adequately  transport  an 
elevated  level  of  modified  turbulent  viscosity  from  the  inlet  boundary  to  the  clutter 
to  model  the  effect  of  free-stream  turbulence  and  influence  separation.  This  aspect  of 
the  Spalart-Allmaras  turbulence  model  is  discussed  in  detail  in  Section  2.3.2 

Because  the  inlet  is  15 D  upstream  of  the  first  row  of  cylinders,  the  high  values  of 
0  set  at  the  inlet  to  mimic  10%  and  15%  turbulence  intensity  (T.I.)  quickly  fall  back 
to  zero  before  reaching  the  clutter.  [4,22]  This  is  clearly  evident  in  Figure  A.l  which 
compares  the  modified  turbulent  viscosity  level  for  the  10%  and  15%  free-stream  T.I. 
cases  along  a  line  describing  y,  z  =  0  down  the  length  of  the  single  cylinder  grid.  The 
cylinder  is  located  at  x  =  0  which  explains  the  absence  of  z>  values  at  x  —  0.  For  the 
15%  T.I.  curve,  the  initial  values  of  z>  are  only  slightly  higher  than  than  their  10%  TI 
counterparts.  In  fact,  velocity  plays  a  much  larger  role  in  elevating  i>.  Regardless, 
prior  to  reaching  the  cylinder,  both  cases  have  been  damped  to  near  zero  levels. 

However,  the  limitation  of  the  Spalart-Allmaras  model  should  not  discourage 
future  researchers  from  attempting  to  assess  the  effect  of  accurately  modeled  free- 
stream  turbulence  intensity.  If  a  DES  is  employed,  then  the  user  should  use  a  k-e  or 
k-u  based  DES  which  would  allow  the  user  to  input  turbulence  intensity  at  the  inlet  via 
direct  specification  of  turbulent  kinetic  energy  and  the  dissipation  rate.  This  option 
is  available  in  the  newest  version  of  Fluent ™:  6.3.21.  Estimation  of  the  turbulent 
kinetic  energy,  k  is  given  by 
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Modified  Turbulent  Viscosity  Along  y  =  0,  z  =  0 
Single  Cylinder  DES 


Figure  A.l:  Modified  Turbulent  Viscosity,  v,  comparison  for  a  single  cylinder  at 

1  m/s. 
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where  I  is  the  turbulence  intensity  given  defined  as 


^ avg 


(A.2) 


By  having  two  equations  describing  the  transport,  generation,  and  destruction 
of  turbulence,  the  user  will  have  more  control  over  the  level  of  modified  turbulent 
viscosity  reaching  the  clutter  element. 


A.2  RANS  to  DES  Transient  Artifacts 

As  noted  in  Section  5.3,  transient  artifacts  persisted  in  the  DES  solution  posing 
as  a  distracting  red  herring.  Figure  ??  depicts  regions  of  positive  and  negative  span- 
wise  velocities  with  magnitudes  upwards  of  50%  of  the  local  stream-wise  velocity  in 
the  cylinder  wake.  This  data  was  collected  after  well  over  a  second  of  DES  which 
hinted  towards  real  flow  physics.  However  after  running  the  simulations  for  three 
more  seconds,  these  periodic  pockets  of  span-wise  flow  vanished  from  the  solution. 
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-0.436  -0.253  -0.0705  0.112  0.295  0.478  0.660 
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-2.30 
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Figure  A. 2:  (a)  Average  ^-velocity  (m/s)  with  5  m/s  (14M  cells). 

x- velocity  (m/s)  with  5  m/s  (14M  cells). 


(b)  Average 
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